#include "stdafx.h"
#include "volumerenderer_private.h"

/*!
\file
This is the main VolumeRenderer implementation file. Initialization related stuff is in volumerenderer_init.cpp.
\sa volumerenderer_private.h
*/


#define VB_NO_GLPROFILE

QSurfaceFormat VolumeRenderer::getDesiredFormat() const
{
	QSurfaceFormat f;

	f.setRenderableType(QSurfaceFormat::OpenGL);
#ifndef VB_NO_GLPROFILE
	//try to get a 3.3 core profile
	f.setProfile(QSurfaceFormat::CoreProfile);
	f.setVersion(3,3);
#endif
	f.setSwapBehavior(QSurfaceFormat::DoubleBuffer);

	//we want 32 bit rgba, 24 bit depth, no stencil is required (probably still provided)
	f.setRedBufferSize(8);
	f.setGreenBufferSize(8);
	f.setBlueBufferSize(8);
	f.setAlphaBufferSize(8);

	//MSAA only applied to boundingbox, lines...
	//volume rendering does not use multisample
	f.setSamples(16);
	f.setDepthBufferSize(24);
	//stencil is -1 by default
	//f.setStencilBufferSize(-1);

#ifdef _DEBUG
	f.setOption(QSurfaceFormat::FormatOption::DebugContext);
#endif

	return f;
}

float VolumeRenderer::getFPS() const
{
	//room for improvement here, current fps calculation method is quite inaccurate
	//better would be averaging over several frames
	if(d.frametime==0)
		return 0.0f;
	else
		return (1000.0f*1000.0f*1000.0f) / d.frametime;
}

int VolumeRenderer::getLastSliceCount() const
{
	return d.sliceCount;
}

int VolumeRenderer::getLastRenderedSliceCount() const
{
	return d.renderedSlices;
}

VolumeMetaData VolumeRenderer::getCurrentVolume() const
{
	return d.currentVolume;
}

bool VolumeRenderer::loadVolume(const QString& metadataFile,VolumeMetaData* returnData)
{
	d.timer.restart();
	//load metadata first
	VolumeMetaData vmd;
	bool ok = d.volLoader->getMetadata(metadataFile,vmd);
	if(!ok)
	{
		qWarning()<<"VolumeLoader returned false while loading metadata, aborting loading";
		return false;
	}
	qDebug()<<"Metadata loaded: "<<vmd;
	qDebug()<<"Metadata loaded in "<<d.timer.restart()<<"ms";


	if(!loadVolume(vmd))
		return false;

	if(returnData)
		*returnData = vmd;

	return true;
}

bool VolumeRenderer::loadVolume(VolumeMetaData& metadata)
{
	d.timer.restart();

	QByteArray data = d.volLoader->readVolume(metadata);
	if(data.isEmpty())
	{
		qWarning()<<"VolumeLoader returned no data, aborting";
		return false;
	}

	qDebug()<<"Volume loaded from file in "<<d.timer.restart()<<"ms";

	d.currentVolume = metadata;
	d.currentVolumeData = data;

	//the c++ side is complete, now we require a gl context to load the volume
	emit requestUpdate();
	emit requestDraw();

	return true;
}

void VolumeRenderer::update()
{
	if(!d.currentVolumeData.isEmpty())
	{
		//new volume data to load
		processVolume(d.currentVolumeData);
		d.currentVolumeData.clear();
	}
}

QVector<float> VolumeRenderer::setupGaussian(int& size)
{
	//size of gaussian kernel:
	int kSize = (6.0f * d.gaussSigma - 1.0f);

	//round up to nearest odd number
	kSize = kSize % 2 == 1 ? kSize : kSize+1;
	// make a minimum of 3x3x3 kernel, max 9^3
	if(kSize<3)
		kSize = 3;
	else if (kSize>9)
		kSize = 9;

	QVector<float> data (kSize*kSize*kSize);

	//inital offsets
	int off_x,off_y,off_z = -(kSize - 1) / 2;

	float denom = 2.0f * (d.gaussSigma * d.gaussSigma);
	float sum = 0.0f;

	//loop da loop
	for(int k = 0; k < kSize; k++) // z slice
	{
		off_y = -(kSize-1)/2;
		for(int j = 0; j < kSize; j++) // y row
		{
			off_x = -(kSize-1)/2;
			for(int i = 0; i < kSize; i++) // x value
			{
				int idx = k*(kSize*kSize) + j*(kSize) + i;

				float num = off_x * off_x + off_y * off_y + off_z * off_z;
				data[idx] = exp(-num / denom);
				sum += data[idx];

				off_x = off_x + 1;
			}
			off_y = off_y + 1;
		}
		off_z = off_z + 1;
	}

	//normalize kernel
	for(int i = 0;i<kSize*kSize*kSize;++i)
	{
		data[i] /= sum;
	}

	size = kSize;
	return data;
}

void VolumeRenderer::processVolume(const QByteArray& data)
{
	qDebug()<<"Initializing volume 3D texture";
	d.timer.restart();

	//remove old volume, if existing
	//(re)create 3d texture with gradients and intensity (used in rendering)
	auto gradTex = d.textures.volume;
	if(gradTex->isCreated())
	{
		gradTex->destroy();	
	}

	//load intensity to 3d texture
	//this texture is only temporary, it is deleted after gradient calculation
	//we have to recreate the 3d texture each time a new volume is loaded, due to qt automatically trying to use immutable storage (this is a good thing)
	QOpenGLTexture volTex(QOpenGLTexture::Target::Target3D);

	//no mipmapping
	volTex.setMipLevelRange(0,0);
	volTex.setMipLevels(0);
	//intensity gets NEAREST filtering
	volTex.setMinMagFilters(QOpenGLTexture::Filter::Nearest,QOpenGLTexture::Filter::Nearest);
	volTex.setSize(d.currentVolume.width,d.currentVolume.height,d.currentVolume.depth);
	volTex.setWrapMode(QOpenGLTexture::WrapMode::ClampToBorder);
	volTex.setBorderColor(0,0,0,0);

	//map alpha to red, others to zero
	volTex.setSwizzleMask(
		QOpenGLTexture::ZeroValue,
		QOpenGLTexture::ZeroValue,
		QOpenGLTexture::ZeroValue,
		QOpenGLTexture::RedValue);

	//set data type and upload
	{
		QOpenGLTexture::PixelType pixType;
		if(d.currentVolume.type == VolumeMetaData::VolumeDataType::ushort)
		{
			volTex.setFormat(QOpenGLTexture::TextureFormat::R16_UNorm);
			pixType = QOpenGLTexture::PixelType::UInt16;
		}
		else
		{
			volTex.setFormat(QOpenGLTexture::TextureFormat::R8_UNorm);
			pixType = QOpenGLTexture::PixelType::UInt8;
		}

		volTex.allocateStorage();

		// set alignment to 1, this is important!
		QOpenGLPixelTransferOptions opts;
		opts.setAlignment(1);
		volTex.setData(QOpenGLTexture::PixelFormat::Red,pixType,(void*)data.data(), &opts);
	}

	volTex.bind(TEXUNIT_INTENSITY);
	qDebug()<<"Intensity volume uploaded to GL in "<<d.timer.restart()<<"ms";


	//*** PREPROCESSING

	QOpenGLShaderProgram* smoothShader = d.shaders.gaussian;
	QOpenGLShaderProgram* gradShader = d.shaders.gradient;

	//use 8x8x8 compute shaders for now, handle uneven volumes later
	GLuint numX = d.currentVolume.width / 8, numY = d.currentVolume.height / 8, numZ = d.currentVolume.depth / 8;

	QOpenGLTexture smoothTex(QOpenGLTexture::Target::Target3D);

	if(d.smoothVolume)
	{
		//create texture for smoothing result
		smoothTex.setMipLevelRange(0,0);
		smoothTex.setMipLevels(0);
		smoothTex.setMinMagFilters(QOpenGLTexture::Filter::Nearest,QOpenGLTexture::Filter::Nearest);
		smoothTex.setSize(d.currentVolume.width,d.currentVolume.height,d.currentVolume.depth);
		smoothTex.setWrapMode(QOpenGLTexture::WrapMode::ClampToBorder);
		smoothTex.setBorderColor(0,0,0,0);

		smoothTex.setSwizzleMask(
			QOpenGLTexture::ZeroValue,
			QOpenGLTexture::ZeroValue,
			QOpenGLTexture::ZeroValue,
			QOpenGLTexture::RedValue);

		//use single float component
		smoothTex.setFormat(QOpenGLTexture::TextureFormat::R32F);
		smoothTex.allocateStorage();
		smoothTex.bind(TEXUNIT_SMOOTH);

		glBindImageTexture(IMGUNIT_SMOOTH,smoothTex.textureId(),0,GL_FALSE,0,GL_WRITE_ONLY,GL_R32F);

		smoothShader->bind();
		smoothShader->setUniformValue("texVolume",TEXUNIT_INTENSITY);
		smoothShader->setUniformValue("imgTarget",IMGUNIT_SMOOTH);

		//calculate gaussian kernel
		int kSize,off;
		QVector<float> data = setupGaussian(kSize);
		qDebug()<<"Using gaussian kernel of size "<<kSize;
		Q_ASSERT(kSize < 10);
		off = - (kSize -1) / 2;
		smoothShader->setUniformValueArray("afWeights",data.constData(),data.size(),1);
		int loc = smoothShader->uniformLocation("ivSize");
		glUniform3i(loc,kSize,kSize,kSize);
		loc = smoothShader->uniformLocation("ivOffset");
		glUniform3i(loc, off,off,off);


		//perform smoothing
		glDispatchCompute(numX,numY,numZ);

		//unbind for cleanliness
		glBindImageTexture(IMGUNIT_SMOOTH,0,0,GL_FALSE,0,GL_WRITE_ONLY,GL_R32F);

		smoothShader->release();
	}

	//setup final volume texture
	gradTex->setMipLevelRange(0,0);
	gradTex->setMipLevels(0);
	//linear filtering here
	gradTex->setMinMagFilters(QOpenGLTexture::Filter::Linear,QOpenGLTexture::Filter::Linear);
	//righthanded transformation in shader, exchange YZ sizes
	gradTex->setSize(d.currentVolume.width,d.currentVolume.depth,d.currentVolume.height);
	//gradTex->setSize(d.currentVolume.width,d.currentVolume.height,d.currentVolume.depth);
	gradTex->setWrapMode(QOpenGLTexture::WrapMode::ClampToEdge);

	//use full floating points for awesomeness
	gradTex->setFormat(QOpenGLTexture::TextureFormat::RGBA32F);
	gradTex->allocateStorage();
	gradTex->bind(TEXUNIT_VOLUME);

	//setup gradTex for image write in gradient shader
	glBindImageTexture(IMGUNIT_VOLUME,gradTex->textureId(),0,GL_FALSE,0,GL_WRITE_ONLY,GL_RGBA32F);

	gradShader->bind();
	if(d.smoothVolume)
		gradShader->setUniformValue("texVolume", TEXUNIT_SMOOTH);
	else
		gradShader->setUniformValue("texVolume", TEXUNIT_INTENSITY);
	gradShader->setUniformValue("imgTarget", IMGUNIT_VOLUME);

	//setup rescaling of values
	uchar wordsize=0u;
	switch(d.currentVolume.type)
	{
	case VolumeMetaData::VolumeDataType::ushort:
		wordsize=16u;
		break;
	case VolumeMetaData::VolumeDataType::uchar:
	default:
		wordsize=8u;
	}
	float scale = 1.0f;
	if(d.currentVolume.actualBits > 0 && d.currentVolume.actualBits<wordsize)
	{
		scale = pow(2,d.currentVolume.actualBits) / pow(2,wordsize);
		qDebug()<<"Data is rescaled to "<<d.currentVolume.actualBits<<"bit range, using "<<scale<<" scale";
		gradShader->setUniformValue("fMaxValue",scale);
	}
	else
	{
		d.currentVolume.actualBits = wordsize;
		gradShader->setUniformValue("fMaxValue",1.0f);
	}

	//fire up the gradient calculation
	glDispatchCompute(numX, numY, numZ);

	gradShader->release();
	//for better reporting of calculation time, we wait until pipeline is finished
	//should only have neglible impact on total performance
	glFinish();
	qDebug()<<"Smoothing and gradient calculation using compute shader in "<<d.timer.restart()<<"ms";

	//unbind for cleanliness
	glBindImageTexture(IMGUNIT_VOLUME,0,0,GL_FALSE,0,GL_WRITE_ONLY,GL_RGBA32F);

	d.volumeLoaded = true;
}

void VolumeRenderer::drawCurrent()
{
	//background is some soft blue color
	//maybe add a gradient in the future
	glClearColor(140/255.0f,143/255.0f,203/255.0f,1.0f);
	//glClearColor(0,0,0,1);
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

	//we dont need no face culling
	//glEnable(GL_CULL_FACE);

	//update matrices
	d.view.setToIdentity();
	d.view.translate(0, 0, -d.distance);
	d.view.rotate(d.rotX, 1, 0, 0);
	d.view.rotate(d.rotY, 0, 1, 0);
	d.mvp = d.mProj * d.view;

	//setup lightvec in world space
	QMatrix4x4 lightRotMat;
	lightRotMat.rotate(d.lightAz,0,1,0);
	lightRotMat.rotate(d.lightAngle,0,0,1);
	d.lightVec = lightRotMat * QVector3D(1,0,0);

	//calc halfway vector
	QVector3D pos = d.view.row(3).toVector3D();
	d.halfVec = (d.lightVec + pos).normalized();

	//draw the current volume
	if(d.volumeLoaded)
	{
		drawVolume();
	}

	//draw optional bounding box
	if(d.renderBB)
	{
		drawBoundingBox();
	}

	if(d.showLightDir)
		drawLine(QVector3D(0,0,0),d.lightVec,QVector3D(1,0,0));

	d.frametime = d.timer.nsecsElapsed();
	d.timer.restart();
}

// Line rendering in this way is hardcore inefficient, only use for debugging
void VolumeRenderer::drawLine(const QVector3D& pointA,const QVector3D& pointB, const QVector3D& color)
{
	//glDisable(GL_DEPTH_TEST);
	glDisable(GL_BLEND);

	d.vaos.line->bind();

	//stuff vertices into buffer
	d.buffers.lineBuffer->bind();
	QVector3D* buf = (QVector3D*)glMapBuffer(GL_ARRAY_BUFFER,GL_WRITE_ONLY);
	buf[0] = pointA;
	buf[1] = pointB;
	glUnmapBuffer(GL_ARRAY_BUFFER);
	d.buffers.lineBuffer->release();

	d.shaders.color->bind();
	d.shaders.color->setUniformValue("mMVP", d.mvp);
	d.shaders.color->setUniformValue("vColor", QVector4D(color,1.0f));

	glDrawArrays(GL_LINES,0,2);
	d.shaders.color->release();
	d.vaos.line->release();
}

void VolumeRenderer::drawBoundingBox()
{
	//bounding box requires multiplicative blending, without depth test
	//glBlendFunc(GL_ZERO, GL_SRC_COLOR);
	glEnable(GL_BLEND);
	glBlendFunc(GL_ZERO,GL_SRC_COLOR);
	//glDisable(GL_DEPTH_TEST);


	d.vaos.box->bind();
	d.shaders.bBox->bind();
	d.shaders.bBox->setUniformValue("mMVP", d.mvp);
	d.shaders.bBox->setUniformValue("colMult", 1.0f);


	glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
	glDrawElements(GL_LINES,2 * Cube::lineCount, GL_UNSIGNED_SHORT, (const void*) d.buffers.lineIndexOffset);

	/*
	//this is for debugging of the bbox shader
	glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
	d.shaders.bBox->setUniformValue("colMult", 0.0f);
	glDrawElements(GL_LINES,2 * Cube::lineCount, GL_UNSIGNED_SHORT, (const void*) d.buffers.lineIndexOffset);
	*/

	d.vaos.box->release();
}

void VolumeRenderer::drawVolume()
{
	//glDisable(GL_DEPTH_TEST);

	d.vaos.slice->bind();
	d.shaders.slice->bind();
	d.sliceCount = 0;
	d.renderedSlices = 0;

	if(! d.fixSlicingMatrix)
	{
		d.mSlicing = d.view;
		d.shaders.slice->setUniformValue("mProj", d.mProj);
	}
	else
	{
		//this is a hack to allow viewing of the slices from all angles when fixing them
		d.shaders.slice->setUniformValue("mProj", d.mProj * d.view * d.mSlicing.inverted());
	}

	//the shader needs the inverted matrix to find the texture coords
	d.shaders.slice->setUniformValue("mInvMV", d.mSlicing.inverted());
	d.shaders.slice->setUniformValue("vWindow",d.window);
	d.shaders.slice->setUniformValue("texVolume",TEXUNIT_VOLUME);
	d.shaders.slice->setUniformValue("vLightVec",d.lightVec);
	d.shaders.slice->setUniformValue("vHalfVec",d.halfVec);

	//transform box into view space
	QVector3D min,max;
	auto transformed = d.bBox.transform3D(d.mSlicing,min,max);

	// the actual range of Z used
	float zRange = max.z() - min.z();
	//total number of slices
	d.sliceCount = qCeil(zRange / d.sliceDelta);

	d.buffers.sliceVBuffer->bind();
	if(!d.doDOF)
	{
		glEnable(GL_BLEND);
		glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

		traverseSlices(transformed,d.sliceCount,min.z(),d.sliceDelta);
	}
	else
	{
		//DOF path

		//the read FB is the default FB, use for initializing with blit
		glBindFramebuffer(GL_READ_FRAMEBUFFER, 0);
		//clearcolor for intermediate buffers is all zero
		glClearColor(0.0f,0.0f,0.0f,0.0f);

		//setup texunits beforehand
		d.shaders.blur->bind();
		d.shaders.blur->setUniformValue("texSlice",TEXUNIT_SLICE);
		d.shaders.blur->setUniformValue("texPrevious",TEXUNIT_PREVIOUS);

		//calc Z of focal plane
		float frontRange = zRange * d.focalPlaneDist;
		float focusZ = max.z() - frontRange;
		float frontDist = max.z() - focusZ;

		int frontSliceCount = qCeil(frontRange / d.sliceDelta);
		int backSliceCount = d.sliceCount - frontSliceCount;

		//traverse back slices BTF

		glBindFramebuffer(GL_DRAW_FRAMEBUFFER,d.fbos.intermediate[1]->handle());
		//fill initial read buffer with current framebuffer contents, this should lead to same blending as normal method
		glBlitFramebuffer(0,0,d.winWidth,d.winHeight,0,0,d.winWidth,d.winHeight,GL_COLOR_BUFFER_BIT,GL_NEAREST);
		glBlendFunc(GL_ONE_MINUS_DST_ALPHA,GL_DST_ALPHA);

		d.shaders.blur->setUniformValue("bFTB",false);
		int intIdx = traverseSlicesDOF(transformed,backSliceCount,min.z(),d.sliceDelta, focusZ, d.fbos.intermediate);
		glActiveTexture(GL_TEXTURE0 + TEXUNIT_INTERMEDIATE);
		glBindTexture(GL_TEXTURE_2D, d.fbos.intermediate[intIdx]->texture());

		//traverse front slices FTB
		glBindFramebuffer(GL_DRAW_FRAMEBUFFER,d.fbos.eye[1]->handle());
		//glClearColor(140/255.0f,143/255.0f,203/255.0f,0.0f);
		glClear(GL_COLOR_BUFFER_BIT);
		//glClearColor(0.0f,0.0f,0.0f,0.0f);
		glBlendFuncSeparate(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA, GL_ONE, GL_ONE_MINUS_SRC_ALPHA);
		d.shaders.blur->bind();
		d.shaders.blur->setUniformValue("bFTB",true);
		int eyeIdx = traverseSlicesDOF(transformed,frontSliceCount,max.z(),-d.sliceDelta, focusZ, d.fbos.eye);

		//blend results together

		//bind the default framebuffer
		glBindFramebuffer(GL_FRAMEBUFFER, 0);

		//fullscreen merge, replacing contents of default FBO
		glDisable(GL_BLEND);
		d.shaders.merge->bind();
		d.shaders.merge->setUniformValue("texIntermediate",TEXUNIT_INTERMEDIATE);
		d.shaders.merge->setUniformValue("texEye",TEXUNIT_PREVIOUS);
		glDrawArrays(GL_TRIANGLES,0,3);		
	}
	d.buffers.sliceVBuffer->release();

	d.vaos.slice->release();
}

/*!
Like \ref traverseSlices, but applies the DOF effects. The returned int is the index of the last draw buffer which should be merged with the other direction result.
*/
int VolumeRenderer::traverseSlicesDOF(const QVector<QVector3D>& transformedBox, int count,float zStart, float zDelta,float focusZ,QOpenGLFramebufferObject** buffers)
{
	GLint offset;
	GLsizei numPoints;

	int drawBuffer = 0;
	int readBuffer = 1;

	//prefetch ids here, qt performs some sync which adds overhead we dont want
	GLuint fbos[2] = {buffers[0]->handle(), buffers[1]->handle()};
	GLuint tex[2] = {buffers[0]->texture(), buffers[1]->texture()};

	GLuint sliceFBO = d.fbos.slice->handle();

	//set and clear the initial read buffer
	glActiveTexture(GL_TEXTURE0 + TEXUNIT_PREVIOUS);
	glBindTexture(GL_TEXTURE_2D,tex[readBuffer]);

	//"reverse" blending factors - basically, destination swapped with source


	for(uint i =0;i<count;++i)
	{
		//move from minimal z (farthest away) to maximal z (nearest)
		float zPos = zStart + i * zDelta;

		//perform slicing and upload slice to GL, if it intersects
		if(updateSliceGeometry(transformedBox,zPos,offset,numPoints))
		{
			++d.renderedSlices;

			//delete old sync object if existed
			if(d.sliceSync[d.currentSliceOffset]!=nullptr)
			{
				glDeleteSync(d.sliceSync[d.currentSliceOffset]);
			}

			//*** SLICE pass

			//bind slice buffer
			//d.fbos.slice->bind();
			//clear slice buffer
			if(zDelta<0)
			{
				//FTB, have to use slice buffer
				glBindFramebuffer(GL_DRAW_FRAMEBUFFER,sliceFBO);
			}
			else
			{
				//BTF, can use intermediate buffer directly with reverse blending
				glBindFramebuffer(GL_DRAW_FRAMEBUFFER,fbos[drawBuffer]);
			}
			//buffers[drawBuffer]->bind();
			glClear(GL_COLOR_BUFFER_BIT);
			glDisable(GL_BLEND);

			//draw the slice
			d.shaders.slice->bind();
			glDrawArrays(GL_TRIANGLE_FAN,offset,numPoints);

			//create sync object for slice
			d.sliceSync[d.currentSliceOffset] = glFenceSync(GL_SYNC_GPU_COMMANDS_COMPLETE,0);

			//*** BLUR pass
			//bind & clear the merge FBO
			if(zDelta<0)
			{
				//blending performed by shader in FTB
				glBindFramebuffer(GL_DRAW_FRAMEBUFFER,fbos[drawBuffer]);
			}
			else
			{
				//blend backwards into the framebuffer
				glEnable(GL_BLEND);
			}
			//buffers[drawBuffer]->bind();
			//glClear(GL_COLOR_BUFFER_BIT);

			//calculate radius of the view-space circle of confusion
			// d = |z_f - z|
			// C(d) = R*d
			float C = d.blurAmount * abs(focusZ - zPos);
			QVector4D cocViewSpace(C,C,zPos,1.0f);
			//transform to clip space
			QVector4D cocClipSpace = d.mProj * cocViewSpace;
			//const QVector2D offset = QVector2D(.5f,.5f);
			QVector2D cocTexSpace = cocClipSpace.toVector2DAffine() * 0.005f;


			//draw FS quad, combining current slice with previous texture
			d.shaders.blur->bind(); 
			d.shaders.blur->setUniformValue("vCoC",cocTexSpace);

			glDrawArrays(GL_TRIANGLES,0,3);

			//bind current tex to PREVIOUS texunit for next iteration
			glBindTexture(GL_TEXTURE_2D,tex[drawBuffer]);		

			//update slice offset
			d.currentSliceOffset = (d.currentSliceOffset + 1) % SLICEBUFFER_SLICECOUNT;

			//swap buffers
			readBuffer = (readBuffer + 1) & 1;
			drawBuffer = (drawBuffer + 1) & 1;
		}
	}

	return readBuffer;
}

/*!
Traverses the transformed bounding box from the specified start point in the specified direction (delta).
*/
void VolumeRenderer::traverseSlices(const QVector<QVector3D>& transformedBox, int count,float zStart, float zDelta)
{
	GLint offset;
	GLsizei numPoints;
	for(uint i =0;i<count;++i)
	{
		//move from minimal z (farthest away) to maximal z (nearest)
		float zPos = zStart + i * zDelta;

		//perform slicing and upload slice to GL, if it intersects
		if(updateSliceGeometry(transformedBox,zPos,offset,numPoints))
		{
			++d.renderedSlices;

			//delete old sync object if existed
			if(d.sliceSync[d.currentSliceOffset]!=nullptr)
			{
				glDeleteSync(d.sliceSync[d.currentSliceOffset]);
			}

			//draw the slice
			glDrawArrays(GL_TRIANGLE_FAN,offset,numPoints);


			//create sync object
			d.sliceSync[d.currentSliceOffset] = glFenceSync(GL_SYNC_GPU_COMMANDS_COMPLETE,0);
			//update slice offset
			d.currentSliceOffset = (d.currentSliceOffset + 1) % SLICEBUFFER_SLICECOUNT;
		}
	}
}

//comparator for std::sort
/*!
A comparator for \c std::sort that uses an external array for lookup and operates on its indices.
*/
template<typename T>
struct ExternalArrayIndexSort
{
	ExternalArrayIndexSort(const T * const data) : data(data)
	{
	}

	bool operator()(int idx1, int idx2)
	{
		//hopefully, you made your data array large enough :P
		return data[idx1]<data[idx2];
	}

	const T * const data;
};


/*!
Performs slicing of a view-aligned plane with the view-transformed volume bounding box, and uploads the slice to OpenGL, if it has intersections.
*/
bool VolumeRenderer::updateSliceGeometry(const QVector<QVector3D>& box, const float zPos, GLint& firstIdx, GLsizei& count)
{
	// Find the intersection points - this is defines our proxy geometry
	// The following is all modeled after http://http.developer.nvidia.com/GPUGems/gpugems_ch39.html, section 39.4.2
	// note that there is a maximum of 6 intersections!

	//setup slicing plane
	const QVector3D planeNormal(0.0f,0.0f,1.0f); //plane is orthogonal to view direction
	QVector3D planePoint(0.0f,0.0f,zPos);

	int numIntersects = 0;
	static QVector3D intersects[8];

	//calculate the centroid on-the-fly to save a loop
	QVector3D center;

	{
		float d;

		//perform X axis intersections
		//these are edges 0-1, 2-3, 4-5, 6-7 of the box
		for(int i =0;i<4;++i)
		{
			const QVector3D& pointA = box.at(i*2);
			const QVector3D& vec = box.at(i*2 + 1) - pointA;

			if(LineIntersectPlane(pointA,vec,planePoint,planeNormal,d))
			{
				intersects[numIntersects] = pointA + d * vec;
				center+=intersects[numIntersects++];
			}
		}

		//perform Y axis intersections
		//these are edges 0-2, 1-3, 4-6, 5-7 of the box
		for(int i =0;i<4;++i)
		{
			const QVector3D& pointA = box.at(2 * i - i % 2); //0 1 4 5...
			const QVector3D& vec = box.at(2 * i - i % 2 + 2) - pointA;
			if(LineIntersectPlane(pointA,vec,planePoint,planeNormal,d))
			{
				intersects[numIntersects] = pointA + d * vec;
				center+=intersects[numIntersects++];
			}
		}

		//perform Z axis intersections
		//these are edges 0-4, 1-5, 2-6, 3-7 of the box
		for(int i =0;i<4;++i)
		{
			const QVector3D& pointA = box.at(i);
			const QVector3D& vec = box.at(i + 4) - pointA;
			if(LineIntersectPlane(pointA,vec,planePoint,planeNormal,d))
			{
				intersects[numIntersects] = pointA + d * vec;
				center+=intersects[numIntersects++];
			}
		}
	}

	//there can not be more than 6 intersections
	Q_ASSERT(numIntersects<=6);	

	//check if there were any usable intersections at all
	if(numIntersects<3)
		return false;

	//finalize center
	center /= numIntersects;

	//find angles around center and init index vector for sorting
	static float angles[6];
	QVector<int> indices(numIntersects); //0,1,2,3.... needed for sorting
	for(int i =0;i<numIntersects;++i)
	{
		indices.replace(i,i);
		const QVector3D& point = intersects[i];
		angles[i] = qAtan2(point.y() - center.y(),point.x() - center.x());
	}

	//perform sorting on the indices
	std::sort(indices.begin(),indices.end(),ExternalArrayIndexSort<float>(angles));

	//update buffer
	//TODO: optimize! use a large buffer with manual synchronization, this should give us faster performance
	GLintptr offset = d.currentSliceOffset * SLICEBUFFER_STRIDE * sizeof(QVector3D);

	//wait until the buffer range is free
	if(d.sliceSync[d.currentSliceOffset]!=nullptr)
	{
		GLenum result = glClientWaitSync(d.sliceSync[d.currentSliceOffset],0,SYNC_TIMEOUT);
		if(result == GL_TIMEOUT_EXPIRED)
		{
			qCritical()<<"GL sync timed out after 1s, aborting slice render";
			return false;
		}
		else if(result == GL_WAIT_FAILED)
		{
			qCritical()<<"GL sync failed, aborting slice render";
			return false;
		}
	}

	//unsynchronized mapping - we perform sync ourselves
	QVector3D* data =(QVector3D*) glMapBufferRange(GL_ARRAY_BUFFER,offset,SLICEBUFFER_STRIDE*sizeof(QVector3D),GL_MAP_WRITE_BIT | GL_MAP_UNSYNCHRONIZED_BIT);
	count = 0;
	if(numIntersects>4) 
	{
		//we only need the center vertex if having more than 4 total vertices, otherwise it would be wasted
		data[count++] = center;
	}
	for(int i =0;i<numIntersects;++i)
	{
		data[count++] = intersects[indices.at(i)];
	}
	if(numIntersects>4)
	{
		//the first vertex must be repeated when in FAN mode
		data[count++] = intersects[indices.at(0)];
	}
	glUnmapBuffer(GL_ARRAY_BUFFER);
	firstIdx = d.currentSliceOffset * SLICEBUFFER_STRIDE;

	return true;
}

void VolumeRenderer::setSize(int width, int height)
{
	d.winWidth = width;
	d.winHeight = height;

	//update viewport
	glViewport(0,0,width,height);

	updateProjection();
	updateFBOs();

}

void VolumeRenderer::setSliceDelta(float sliceDelta)
{
	d.sliceDelta = sliceDelta;
	emit requestDraw();
}

void VolumeRenderer::setLowerWindowLimit(float w)
{
	d.window.setX(w);
	emit requestDraw();
}

void VolumeRenderer::setUpperWindowLimit(float w)
{
	d.window.setY(w);
	emit requestDraw();
}

void VolumeRenderer::setProjectionType(VolumeRenderer::ProjectionType pt)
{
	d.projType = pt;
	updateProjection();

	emit requestDraw();
}

void VolumeRenderer::updateProjection()
{
	d.mProj.setToIdentity();
	switch (d.projType)
	{
	case ProjectionType::Orthographic:
		d.mProj.ortho(-1,1,-1,1,0,10);
		break;
	case ProjectionType::Perspective:
	default:
		d.mProj.perspective(60,(float)d.winWidth/d.winHeight,0.1,100.0);
		break;
	}

}

void VolumeRenderer::setLightParameters(float azimuth, float angle)
{
	d.lightAz = azimuth;
	d.lightAngle = angle;

	const float xExtent = 90.0f;

	//clamp angle to [-90,90]
	if(d.lightAngle>xExtent)
		d.lightAngle=xExtent;
	if(d.lightAngle<-xExtent)
		d.lightAngle=-xExtent;

	//normalize y rotation to [0,360]
	while(d.lightAz>360.0f)
		d.lightAz-=360.0f;
	while(d.lightAz<0.0f)
		d.lightAz+=360.0f;

	emit lightParamsChanged(d.lightAz,d.lightAngle);
	emit requestDraw();
}

void VolumeRenderer::lightDelta(float azimuth, float angle)
{
	setLightParameters(d.lightAz + azimuth, d.lightAngle + angle);
}

void VolumeRenderer::setRotation(float x,float y)
{
	d.rotX = x;
	d.rotY = y;
	//const float xExtent = 90.0f;

	//clamp x rotation
	//if(d.rotX>xExtent)
	//	d.rotX=xExtent;
	//if(d.rotX<-xExtent)
	//	d.rotX=-xExtent;

	//normalize x rot
	while(d.rotX>360.0f)
		d.rotX-=360.0f;
	while(d.rotX<0.0f)
		d.rotX+=360.0f;

	//normalize y rotation
	while(d.rotY>360.0f)
		d.rotY-=360.0f;
	while(d.rotY<0.0f)
		d.rotY+=360.0f;

	emit rotationChanged(d.rotX, d.rotY);
	emit requestDraw();
}

void VolumeRenderer::rotationDelta(float deltaX,float deltaY)
{
	setRotation(d.rotX+deltaX,d.rotY+deltaY);
}

void VolumeRenderer::setDistance(float dst)
{
	d.distance = dst;

	const float maxDist=10.0f, minDist = 0.5f;
	if(d.distance > maxDist)
		d.distance = maxDist;
	if(d.distance < minDist)
		d.distance = minDist;

	emit distanceChanged(d.distance);
	emit requestDraw();
}

void VolumeRenderer::distanceDelta(float dst)
{
	setDistance(d.distance + dst);
}

void VolumeRenderer::setRenderBoundingBox(bool b)
{
	d.renderBB = b;
	emit requestDraw();
}

void VolumeRenderer::setFixSlicingMatrix(bool b)
{
	d.fixSlicingMatrix = b;
	emit requestDraw();
}

void VolumeRenderer::setShowLightDirection(bool v)
{
	d.showLightDir = v;
	emit requestDraw();
}

void VolumeRenderer::setSmoothVolume(bool v)
{
	d.smoothVolume = v;

	//reload volume if one is loaded
	if(d.volumeLoaded)
	{
		qDebug()<<"Smoothing change requires reload...";
		VolumeMetaData vmd = d.currentVolume;
		loadVolume(vmd);
	}
}

void VolumeRenderer::setGaussianSigma(float s)
{
	d.gaussSigma = s;
	//reload volume if one is loaded
	if(d.volumeLoaded)
	{
		qDebug()<<"Smoothing change requires reload...";
		VolumeMetaData vmd = d.currentVolume;
		loadVolume(vmd);
	}
}

void VolumeRenderer::setApplyDOF(bool v)
{
	d.doDOF = v;
	emit requestDraw();
}

void VolumeRenderer::setBlurAmount(float val)
{
	d.blurAmount = val;
	emit requestDraw();
}

void VolumeRenderer::setFocalPlane(float val)
{
	d.focalPlaneDist=val;
	emit requestDraw();
}