
#include "itkCUDARayCast.cuh"

#include <cuda.h>
#include <cuda_runtime_api.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "texture_fetch_functions.h"
#include <algorithm> 
#include <iostream>
#include <vector>

#define KERNEL                      __global__
#define HOST                        __host__
#define DEVICE                      __device__
#define HOST_AND_DEVICE             __host__ __device__
#define DEVICE_CONST                __device__ __constant__

// This variable contains the DICOM set
float* d_object3D;

// This variable contains the 2D output from CUDA
float *d_object2D;

// Constants depending on the DICOM
DEVICE_CONST int d_sizeCT[3];
DEVICE_CONST  float ctPixelSpacing[3];

// Constant depending on image output
DEVICE_CONST int DRRImageSize[2];

// Constants dependion on the specific DRR
DEVICE_CONST  float d_DRR_Parameters[12];

// This variable contains the DICOM loaded as a Texture ( read-only, fast-cached memory)
cudaTextureObject_t tex_object3D = 0;

cudaStream_t stream1;


void loadDICOMInGPUMemory(float *cpp_object3D, int *sizeCT, float *pixelSpacing)
{
	long int object3Dsize = sizeCT[0] * sizeCT[1] * sizeCT[2];
	cudaStreamCreateWithFlags(&stream1, cudaStreamNonBlocking);


	cudaMalloc((void**)&d_object3D, object3Dsize * sizeof(float));
	cudaMemcpyAsync(d_object3D, cpp_object3D, object3Dsize * sizeof(float), cudaMemcpyHostToDevice, stream1);

	cudaMemcpyToSymbol(ctPixelSpacing, pixelSpacing, 3 * sizeof(float), 0, cudaMemcpyHostToDevice);
	cudaMemcpyToSymbol(d_sizeCT, sizeCT, 3 * sizeof(int), 0, cudaMemcpyHostToDevice);


	// create texture object
	cudaResourceDesc resDesc;
	memset(&resDesc, 0, sizeof(resDesc));
	resDesc.resType = cudaResourceTypeLinear;
	resDesc.res.linear.devPtr = d_object3D;
	resDesc.res.linear.desc.f = cudaChannelFormatKindFloat;
	resDesc.res.linear.desc.x = 32; // bits per channel
	resDesc.res.linear.sizeInBytes = object3Dsize * sizeof(float);

	cudaTextureDesc texDesc;
	memset(&texDesc, 0, sizeof(texDesc));
	texDesc.readMode = cudaReadModeElementType;

	// create CUDA texture object
	cudaDestroyTextureObject(tex_object3D);
	cudaCreateTextureObject(&tex_object3D, &resDesc, &texDesc, NULL);

	cudaStreamDestroy(stream1);

}


void loadOuputVariablesInGPUMemory(int dimX, int dimZ)
{
	long int vectorSize = dimX * dimZ;
	int OutputImageSize[2] = { dimX, dimZ };

	cudaMalloc((void**)&d_object2D, vectorSize * sizeof(float));

	cudaMemcpyToSymbol(DRRImageSize, OutputImageSize, 2 * sizeof(int), 0, cudaMemcpyHostToDevice);
}


void freeDICOMFromGPUMemory()
{
	cudaFree(d_object3D);
}

void freeAuxiliaryVariablesInGPUMemory()
{
	cudaFree(d_object2D);
}


__global__ void drrCUDA(float* object2D, cudaTextureObject_t tex_object3D)
{

	float stepInX[3] = { d_DRR_Parameters[0],d_DRR_Parameters[1],d_DRR_Parameters[2] };
	float stepInY[3] = { d_DRR_Parameters[3],d_DRR_Parameters[4],d_DRR_Parameters[5] };
	float corner00[3] = { d_DRR_Parameters[6],d_DRR_Parameters[7],d_DRR_Parameters[8] };
	float SourceWorld[3] = { d_DRR_Parameters[9],d_DRR_Parameters[10],d_DRR_Parameters[11] };

	int total_dx = DRRImageSize[0];
	int total_dz = DRRImageSize[1];

	//Every thread calculates its own id number
	long int idx = (blockIdx.x*blockDim.x) + threadIdx.x;

	// This checks if the thread number is bigger than the amount of pixels
	if (idx >= total_dx * total_dz)
		return;

	// Converting number of pixels to rows and columns
	int dz = idx / total_dx;
	int dx = idx - dz * total_dx;

	//Calculate the spatial position of the pixel
	float drrPixelWorld[3] = { 0 };
	drrPixelWorld[0] = corner00[0] + ((stepInX[0])*dx) + ((stepInY[0])*dz);
	drrPixelWorld[1] = corner00[1] + ((stepInX[1])*dx) + ((stepInY[1])*dz);
	drrPixelWorld[2] = corner00[2] + ((stepInX[2])*dx) + ((stepInY[2])*dz);

	//Calculate the ray vector
	float rayVector[3] = { 0 };
	rayVector[0] = drrPixelWorld[0] - SourceWorld[0];
	rayVector[1] = drrPixelWorld[1] - SourceWorld[1];
	rayVector[2] = drrPixelWorld[2] - SourceWorld[2];

	double pixval = 0;
	// need
	double transformedFocalPoint[3];
	double rayDirection[3];
	double origin[3];
	double spacing[3];
	double size[3];
	double rayPosition[3];
	double threshold;

	/////////
	bool m_ValidRay = false;
	double m_RayVoxelStartPosition[3];
	double m_RayVoxelEndPosition[3];
	double m_Position3Dvox[3];
	double m_VoxelIncrement[3];
	int m_TraversalDirection;
	int m_TotalRayVoxelPlanes;
	int m_NumVoxelPlanesTraversed;
	double* m_RayIntersectionVoxels[4];
	int m_RayIntersectionVoxelIndex[3];
	int m_NumberOfVoxelsInx = 0;
	int m_NumberOfVoxelsIny = 0;
	int m_NumberOfVoxelsInz = 0;
	double m_VoxelDimensionInx = 0;
	double m_VoxelDimensionIny = 0;
	double m_VoxelDimensionInz = 0;
	double m_RayStartCoordInMM[3];
	double m_RayEndCoordInMM[3];
	double m_BoundingPlane[6][4];
	double m_BoundingCorner[8][3];
	double m_CurrentRayPositionInMM[3];
	double m_RayDirectionInMM[3];

	//zero state
	for (int i = 0; i < 3; i++) {
		m_CurrentRayPositionInMM[i] = 0.;
		m_RayDirectionInMM[i] = 0.;
		m_RayVoxelStartPosition[i] = 0.;
		m_RayVoxelEndPosition[i] = 0.;
		m_VoxelIncrement[i] = 0.;
		m_RayIntersectionVoxelIndex[i] = 0;
	}
	m_TraversalDirection = 0;
	m_TotalRayVoxelPlanes = 0;
	m_NumVoxelPlanesTraversed = -1;
	for (int i = 0; i < 4; i++) {
		m_RayIntersectionVoxels[i] = nullptr;
	}

	//initialise
	/*
		recordVolumeDimension
		DefineCorners
		CalcPlanesAndCorners
	*/
	m_NumberOfVoxelsInx = size[0];
	m_NumberOfVoxelsIny = size[1];
	m_NumberOfVoxelsInz = size[2];
	m_VoxelDimensionInx = spacing[0];
	m_VoxelDimensionIny = spacing[1];
	m_VoxelDimensionInz = spacing[2];

	m_BoundingCorner[0][0] = m_BoundingCorner[1][0] = m_BoundingCorner[2][0] =
		m_BoundingCorner[3][0] = 0;
	m_BoundingCorner[4][0] = m_BoundingCorner[5][0] = m_BoundingCorner[6][0] =
		m_BoundingCorner[7][0] = m_VoxelDimensionInx * m_NumberOfVoxelsInx;

	m_BoundingCorner[1][1] = m_BoundingCorner[3][1] = m_BoundingCorner[5][1] =
		m_BoundingCorner[7][1] = m_VoxelDimensionIny * m_NumberOfVoxelsIny;

	m_BoundingCorner[0][1] = m_BoundingCorner[2][1] = m_BoundingCorner[4][1] =
		m_BoundingCorner[6][1] = 0;

	m_BoundingCorner[0][2] = m_BoundingCorner[1][2] = m_BoundingCorner[4][2] =
		m_BoundingCorner[5][2] = m_VoxelDimensionInz * m_NumberOfVoxelsInz;

	m_BoundingCorner[2][2] = m_BoundingCorner[3][2] = m_BoundingCorner[6][2] =
		m_BoundingCorner[7][2] = 0;

	int j;

	// find the equations of the planes

	int c1 = 0, c2 = 0, c3 = 0;

	for (j = 0; j < 6; j++)
	{                                // loop around for planes
		switch (j)
		{                // which corners to take
		case 0:
			c1 = 1; c2 = 2; c3 = 3;
			break;
		case 1:
			c1 = 4; c2 = 5; c3 = 6;
			break;
		case 2:
			c1 = 5; c2 = 3; c3 = 7;
			break;
		case 3:
			c1 = 2; c2 = 4; c3 = 6;
			break;
		case 4:
			c1 = 1; c2 = 5; c3 = 0;
			break;
		case 5:
			c1 = 3; c2 = 7; c3 = 2;
			break;
		}

		double line1x, line1y, line1z;
		double line2x, line2y, line2z;

		// lines from one corner to another in x,y,z dirns
		line1x = m_BoundingCorner[c1][0] - m_BoundingCorner[c2][0];
		line2x = m_BoundingCorner[c1][0] - m_BoundingCorner[c3][0];

		line1y = m_BoundingCorner[c1][1] - m_BoundingCorner[c2][1];
		line2y = m_BoundingCorner[c1][1] - m_BoundingCorner[c3][1];

		line1z = m_BoundingCorner[c1][2] - m_BoundingCorner[c2][2];
		line2z = m_BoundingCorner[c1][2] - m_BoundingCorner[c3][2];

		double A, B, C, D;

		// take cross product
		A = line1y * line2z - line2y * line1z;
		B = line2x * line1z - line1x * line2z;
		C = line1x * line2y - line2x * line1y;

		// find constant
		D = -(A * m_BoundingCorner[c1][0]
			+ B * m_BoundingCorner[c1][1]
			+ C * m_BoundingCorner[c1][2]);

		// initialise plane value and normalise
		m_BoundingPlane[j][0] = A / std::sqrt(A * A + B * B + C * C);
		m_BoundingPlane[j][1] = B / std::sqrt(A * A + B * B + C * C);
		m_BoundingPlane[j][2] = C / std::sqrt(A * A + B * B + C * C);
		m_BoundingPlane[j][3] = D / std::sqrt(A * A + B * B + C * C);
		// A*A+B*B+C*C == 0 »á³ö´í

	}

	//setray
	m_NumberOfVoxelsInx = size[0];
	m_NumberOfVoxelsIny = size[1];
	m_NumberOfVoxelsInz = size[2];

	m_VoxelDimensionInx = spacing[0];
	m_VoxelDimensionIny = spacing[1];
	m_VoxelDimensionInz = spacing[2];

	m_CurrentRayPositionInMM[0] = rayPosition[0];
	m_CurrentRayPositionInMM[1] = rayPosition[1];
	m_CurrentRayPositionInMM[2] = rayPosition[2];

	m_RayDirectionInMM[0] = rayDirection[0];
	m_RayDirectionInMM[1] = rayDirection[1];
	m_RayDirectionInMM[2] = rayDirection[2];

	bool   noInterceptFlag[6];
	double cubeIntercepts[6][3];
	constexpr unsigned int numSides = 6;
	// Calculate intercept of ray with planes
	double interceptx[6];
	double intercepty[6];
	double interceptz[6];
	for (unsigned int j = 0; j < numSides; ++j)
	{
		const double denom = (m_BoundingPlane[j][0] * m_RayDirectionInMM[0]
			+ m_BoundingPlane[j][1] * m_RayDirectionInMM[1]
			+ m_BoundingPlane[j][2] * m_RayDirectionInMM[2]);

		if ((long)(denom * 100) != 0)
		{
			const double d = -(m_BoundingPlane[j][3]
				+ m_BoundingPlane[j][0] * m_CurrentRayPositionInMM[0]
				+ m_BoundingPlane[j][1] * m_CurrentRayPositionInMM[1]
				+ m_BoundingPlane[j][2] * m_CurrentRayPositionInMM[2]) / denom;

			interceptx[j] = m_CurrentRayPositionInMM[0] + d * m_RayDirectionInMM[0];
			intercepty[j] = m_CurrentRayPositionInMM[1] + d * m_RayDirectionInMM[1];
			interceptz[j] = m_CurrentRayPositionInMM[2] + d * m_RayDirectionInMM[2];

			noInterceptFlag[j] = true;  //OK
		}
		else
		{
			noInterceptFlag[j] = false;  //NOT OK
		}
	}
	unsigned int nSidesCrossed = 0;
	for (unsigned int j = 0; j < numSides; ++j)
	{
		// Work out which corners to use

		int c[4];
		if (j == 0)
		{
			c[0] = 0; c[1] = 1; c[2] = 3; c[3] = 2;
		}
		else if (j == 1)
		{
			c[0] = 4; c[1] = 5; c[2] = 7; c[3] = 6;
		}
		else if (j == 2)
		{
			c[0] = 1; c[1] = 5; c[2] = 7; c[3] = 3;
		}
		else if (j == 3)
		{
			c[0] = 0; c[1] = 2; c[2] = 6; c[3] = 4;
		}
		else if (j == 4)
		{ //TOP
			c[0] = 0; c[1] = 1; c[2] = 5; c[3] = 4;
		}
		else if (j == 5)
		{ //BOTTOM
			c[0] = 2; c[1] = 3; c[2] = 7; c[3] = 6;
		}
		// Calculate vectors from corner of ct volume to intercept.
		double cornerVect[4][3];
		for (unsigned int i = 0; i < 4; ++i)
		{
			if (noInterceptFlag[j])
			{
				cornerVect[i][0] = m_BoundingCorner[c[i]][0] - interceptx[j];
				cornerVect[i][1] = m_BoundingCorner[c[i]][1] - intercepty[j];
				cornerVect[i][2] = m_BoundingCorner[c[i]][2] - interceptz[j];
			}
			else
			{
				cornerVect[i][0] = 0;
				cornerVect[i][1] = 0;
				cornerVect[i][2] = 0;
			}
		}
		// Do cross product with these vectors
		int cross[4][3];
		unsigned int k = 0;
		for (unsigned int i = 0; i < 4; ++i)
		{
			if (i == 3)
			{
				k = 0;
			}
			else
			{
				k = i + 1;
			}
			const double ax = cornerVect[i][0];
			const double ay = cornerVect[i][1];
			const double az = cornerVect[i][2];
			const double bx = cornerVect[k][0];
			const double by = cornerVect[k][1];
			const double bz = cornerVect[k][2];

			// The int and divide by 100 are to avoid rounding errors.  If
			// these are not included then you get values fluctuating around
			// zero and so in the subsequent check, all the values are not
			// above or below zero.  NB. If you "INT" by too much here though
			// you can get problems in the corners of your volume when rays
			// are allowed to go through more than one plane.
			cross[i][0] = static_cast<int>((ay * bz - az * by) / 100);
			cross[i][1] = static_cast<int>((az * bx - ax * bz) / 100);
			cross[i][2] = static_cast<int>((ax * by - ay * bx) / 100);
		}

		// See if a sign change occurred between all these cross products
		// if not, then the ray went through this plane

		unsigned int crossFlag = 0;
		for (unsigned int i = 0; i < 3; ++i)
		{
			if ((cross[0][i] <= 0
				&& cross[1][i] <= 0
				&& cross[2][i] <= 0
				&& cross[3][i] <= 0)

				|| (cross[0][i] >= 0
					&& cross[1][i] >= 0
					&& cross[2][i] >= 0
					&& cross[3][i] >= 0))
			{
				++crossFlag;
			}
		}

		if (crossFlag == 3 && noInterceptFlag[j] == 1)
		{
			cubeIntercepts[nSidesCrossed][0] = interceptx[j];
			cubeIntercepts[nSidesCrossed][1] = intercepty[j];
			cubeIntercepts[nSidesCrossed][2] = interceptz[j];
			++nSidesCrossed;
		}
	} // End of loop over all four planes

	m_RayStartCoordInMM[0] = cubeIntercepts[0][0];
	m_RayStartCoordInMM[1] = cubeIntercepts[0][1];
	m_RayStartCoordInMM[2] = cubeIntercepts[0][2];

	m_RayEndCoordInMM[0] = cubeIntercepts[1][0];
	m_RayEndCoordInMM[1] = cubeIntercepts[1][1];
	m_RayEndCoordInMM[2] = cubeIntercepts[1][2];

	if (nSidesCrossed >= 3)
	{
		double maxInterDist = 0.0;
		for (unsigned int j = 0; j < nSidesCrossed - 1; ++j)
		{
			for (unsigned int k = j + 1; k < nSidesCrossed; ++k)
			{
				double interDist = 0.0;
				for (unsigned int i = 0; i < 3; ++i)
				{
					interDist += (cubeIntercepts[j][i] - cubeIntercepts[k][i])
						* (cubeIntercepts[j][i] - cubeIntercepts[k][i]);
				}
				if (interDist > maxInterDist)
				{
					maxInterDist = interDist;

					m_RayStartCoordInMM[0] = cubeIntercepts[j][0];
					m_RayStartCoordInMM[1] = cubeIntercepts[j][1];
					m_RayStartCoordInMM[2] = cubeIntercepts[j][2];

					m_RayEndCoordInMM[0] = cubeIntercepts[k][0];
					m_RayEndCoordInMM[1] = cubeIntercepts[k][1];
					m_RayEndCoordInMM[2] = cubeIntercepts[k][2];
				}
			}
		}
		nSidesCrossed = 2;
	}

	if (nSidesCrossed == 2) {
		m_ValidRay = true;
	}
	else {
		m_ValidRay = false;
	}
	if (!m_ValidRay) {
		m_NumVoxelPlanesTraversed = -1;
		for (int i = 0; i < 3; i++)
		{
			m_RayVoxelStartPosition[i] = 0.;
			m_RayVoxelEndPosition[i] = 0.;
			m_VoxelIncrement[i] = 0.;
			m_RayIntersectionVoxelIndex[i] = 0;
		}
		m_TraversalDirection = 0;

		m_TotalRayVoxelPlanes = 0;

		for (int i = 0; i < 4; i++)
		{
			m_RayIntersectionVoxels[i] = nullptr;
		}
		return;
	}
	//endpointsinvoxels
	m_RayVoxelStartPosition[0] = m_RayStartCoordInMM[0] / m_VoxelDimensionInx;
	m_RayVoxelStartPosition[1] = m_RayStartCoordInMM[1] / m_VoxelDimensionIny;
	m_RayVoxelStartPosition[2] = m_RayStartCoordInMM[2] / m_VoxelDimensionInz;

	m_RayVoxelEndPosition[0] = m_RayEndCoordInMM[0] / m_VoxelDimensionInx;
	m_RayVoxelEndPosition[1] = m_RayEndCoordInMM[1] / m_VoxelDimensionIny;
	m_RayVoxelEndPosition[2] = m_RayEndCoordInMM[2] / m_VoxelDimensionInz;
	//calcdirnvector
	double xNum, yNum, zNum;

	// Calculate the number of voxels in each direction

	xNum = std::fabs(m_RayVoxelStartPosition[0] - m_RayVoxelEndPosition[0]);
	yNum = std::fabs(m_RayVoxelStartPosition[1] - m_RayVoxelEndPosition[1]);
	zNum = std::fabs(m_RayVoxelStartPosition[2] - m_RayVoxelEndPosition[2]);

	// The direction iterated in is that with the greatest number of voxels
	// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	// Iterate in X direction

	if ((xNum >= yNum) && (xNum >= zNum))
	{
		if (m_RayVoxelStartPosition[0] < m_RayVoxelEndPosition[0])
		{
			m_VoxelIncrement[0] = 1;

			m_VoxelIncrement[1] =
				(m_RayVoxelStartPosition[1]
					- m_RayVoxelEndPosition[1]) / (m_RayVoxelStartPosition[0]
						- m_RayVoxelEndPosition[0]);

			m_VoxelIncrement[2] =
				(m_RayVoxelStartPosition[2]
					- m_RayVoxelEndPosition[2]) / (m_RayVoxelStartPosition[0]
						- m_RayVoxelEndPosition[0]);
		}
		else
		{
			m_VoxelIncrement[0] = -1;

			m_VoxelIncrement[1] =
				-(m_RayVoxelStartPosition[1]
					- m_RayVoxelEndPosition[1]) / (m_RayVoxelStartPosition[0]
						- m_RayVoxelEndPosition[0]);

			m_VoxelIncrement[2] =
				-(m_RayVoxelStartPosition[2]
					- m_RayVoxelEndPosition[2]) / (m_RayVoxelStartPosition[0]
						- m_RayVoxelEndPosition[0]);
		}

		// This section is to alter the start position in order to
		// place the center of the voxels in there correct positions,
		// rather than placing them at the corner of voxels which is
		// what happens if this is not carried out.  The reason why
		// x has no -0.5 is because this is the direction we are going
		// to iterate in and therefore we wish to go from center to
		// center rather than finding the surrounding voxels.

		m_RayVoxelStartPosition[1] += ((int)m_RayVoxelStartPosition[0]
			- m_RayVoxelStartPosition[0]) * m_VoxelIncrement[1] * m_VoxelIncrement[0]
			+ 0.5 * m_VoxelIncrement[1] - 0.5;

		m_RayVoxelStartPosition[2] += ((int)m_RayVoxelStartPosition[0]
			- m_RayVoxelStartPosition[0]) * m_VoxelIncrement[2] * m_VoxelIncrement[0]
			+ 0.5 * m_VoxelIncrement[2] - 0.5;

		m_RayVoxelStartPosition[0] = (int)m_RayVoxelStartPosition[0] + 0.5 * m_VoxelIncrement[0];

		m_TotalRayVoxelPlanes = (int)xNum;

		m_TraversalDirection = 1;
	}

	// Iterate in Y direction

	else if ((yNum >= xNum) && (yNum >= zNum))
	{
		if (m_RayVoxelStartPosition[1] < m_RayVoxelEndPosition[1])
		{
			m_VoxelIncrement[1] = 1;

			m_VoxelIncrement[0] =
				(m_RayVoxelStartPosition[0]
					- m_RayVoxelEndPosition[0]) / (m_RayVoxelStartPosition[1]
						- m_RayVoxelEndPosition[1]);

			m_VoxelIncrement[2] =
				(m_RayVoxelStartPosition[2]
					- m_RayVoxelEndPosition[2]) / (m_RayVoxelStartPosition[1]
						- m_RayVoxelEndPosition[1]);
		}
		else
		{
			m_VoxelIncrement[1] = -1;

			m_VoxelIncrement[0] =
				-(m_RayVoxelStartPosition[0]
					- m_RayVoxelEndPosition[0]) / (m_RayVoxelStartPosition[1]
						- m_RayVoxelEndPosition[1]);

			m_VoxelIncrement[2] =
				-(m_RayVoxelStartPosition[2]
					- m_RayVoxelEndPosition[2]) / (m_RayVoxelStartPosition[1]
						- m_RayVoxelEndPosition[1]);
		}

		m_RayVoxelStartPosition[0] += ((int)m_RayVoxelStartPosition[1]
			- m_RayVoxelStartPosition[1]) * m_VoxelIncrement[0] * m_VoxelIncrement[1]
			+ 0.5 * m_VoxelIncrement[0] - 0.5;

		m_RayVoxelStartPosition[2] += ((int)m_RayVoxelStartPosition[1]
			- m_RayVoxelStartPosition[1]) * m_VoxelIncrement[2] * m_VoxelIncrement[1]
			+ 0.5 * m_VoxelIncrement[2] - 0.5;

		m_RayVoxelStartPosition[1] = (int)m_RayVoxelStartPosition[1] + 0.5 * m_VoxelIncrement[1];

		m_TotalRayVoxelPlanes = (int)yNum;

		m_TraversalDirection = 2;
	}

	// Iterate in Z direction

	else
	{
		if (m_RayVoxelStartPosition[2] < m_RayVoxelEndPosition[2])
		{
			m_VoxelIncrement[2] = 1;

			m_VoxelIncrement[0] =
				(m_RayVoxelStartPosition[0]
					- m_RayVoxelEndPosition[0]) / (m_RayVoxelStartPosition[2]
						- m_RayVoxelEndPosition[2]);

			m_VoxelIncrement[1] =
				(m_RayVoxelStartPosition[1]
					- m_RayVoxelEndPosition[1]) / (m_RayVoxelStartPosition[2]
						- m_RayVoxelEndPosition[2]);
		}
		else
		{
			m_VoxelIncrement[2] = -1;

			m_VoxelIncrement[0] =
				-(m_RayVoxelStartPosition[0]
					- m_RayVoxelEndPosition[0]) / (m_RayVoxelStartPosition[2]
						- m_RayVoxelEndPosition[2]);

			m_VoxelIncrement[1] =
				-(m_RayVoxelStartPosition[1]
					- m_RayVoxelEndPosition[1]) / (m_RayVoxelStartPosition[2]
						- m_RayVoxelEndPosition[2]);
		}

		m_RayVoxelStartPosition[0] += ((int)m_RayVoxelStartPosition[2]
			- m_RayVoxelStartPosition[2]) * m_VoxelIncrement[0] * m_VoxelIncrement[2]
			+ 0.5 * m_VoxelIncrement[0] - 0.5;

		m_RayVoxelStartPosition[1] += ((int)m_RayVoxelStartPosition[2]
			- m_RayVoxelStartPosition[2]) * m_VoxelIncrement[1] * m_VoxelIncrement[2]
			+ 0.5 * m_VoxelIncrement[1] - 0.5;

		m_RayVoxelStartPosition[2] = (int)m_RayVoxelStartPosition[2] + 0.5 * m_VoxelIncrement[2];

		m_TotalRayVoxelPlanes = (int)zNum;

		m_TraversalDirection = 3;
	}
	//adjustraylength
	bool startOK, endOK;

	int Istart[3];
	int Idirn[3];

	if (m_TraversalDirection == 1)
	{
		Idirn[0] = 0;
		Idirn[1] = 1;
		Idirn[2] = 1;
	}
	else if (m_TraversalDirection == 2)
	{
		Idirn[0] = 1;
		Idirn[1] = 0;
		Idirn[2] = 1;
	}
	else if (m_TraversalDirection == 3)
	{
		Idirn[0] = 1;
		Idirn[1] = 1;
		Idirn[2] = 0;
	}

	do
	{
		startOK = false;
		endOK = false;

		Istart[0] = (int)std::floor(m_RayVoxelStartPosition[0]);
		Istart[1] = (int)std::floor(m_RayVoxelStartPosition[1]);
		Istart[2] = (int)std::floor(m_RayVoxelStartPosition[2]);

		if ((Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInx)
			&& (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsIny)
			&& (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInz))
		{
			startOK = true;
		}
		else
		{
			m_RayVoxelStartPosition[0] += m_VoxelIncrement[0];
			m_RayVoxelStartPosition[1] += m_VoxelIncrement[1];
			m_RayVoxelStartPosition[2] += m_VoxelIncrement[2];

			m_TotalRayVoxelPlanes--;
		}

		Istart[0] = (int)std::floor(m_RayVoxelStartPosition[0]
			+ m_TotalRayVoxelPlanes * m_VoxelIncrement[0]);

		Istart[1] = (int)std::floor(m_RayVoxelStartPosition[1]
			+ m_TotalRayVoxelPlanes * m_VoxelIncrement[1]);

		Istart[2] = (int)std::floor(m_RayVoxelStartPosition[2]
			+ m_TotalRayVoxelPlanes * m_VoxelIncrement[2]);

		if ((Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInx)
			&& (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsIny)
			&& (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInz))
		{
			endOK = true;
		}
		else
		{
			m_TotalRayVoxelPlanes--;
		}
	} while ((!(startOK && endOK)) && (m_TotalRayVoxelPlanes > 1));

	m_ValidRay = (startOK && endOK);
	//reset
	m_NumVoxelPlanesTraversed = -1;
	if (!m_ValidRay) {
		for (int i = 0; i < 3; i++)
		{
			m_RayVoxelStartPosition[i] = 0.;
			m_RayVoxelEndPosition[i] = 0.;
			m_VoxelIncrement[i] = 0.;
			m_RayIntersectionVoxelIndex[i] = 0;
		}
		m_TraversalDirection = 0;

		m_TotalRayVoxelPlanes = 0;

		for (int i = 0; i < 4; i++)
		{
			m_RayIntersectionVoxels[i] = nullptr;
		}
	}
	else {
		for (int i = 0; i < 3; i++)
		{
			m_Position3Dvox[i] = m_RayVoxelStartPosition[i];
		}

		int Ix, Iy, Iz;

		Ix = (int)(m_RayVoxelStartPosition[0]);
		Iy = (int)(m_RayVoxelStartPosition[1]);
		Iz = (int)(m_RayVoxelStartPosition[2]);

		m_RayIntersectionVoxelIndex[0] = Ix;
		m_RayIntersectionVoxelIndex[1] = Iy;
		m_RayIntersectionVoxelIndex[2] = Iz;

		switch (m_TraversalDirection)
		{
		case 1:
		{
			if ((Ix >= 0) && (Ix < m_NumberOfVoxelsInx)
				&& (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsIny)
				&& (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInz))
			{
				index[0] = Ix; index[1] = Iy; index[2] = Iz;
				m_RayIntersectionVoxels[0] =
					this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index);

				index[0] = Ix; index[1] = Iy + 1; index[2] = Iz;
				m_RayIntersectionVoxels[1] =
					(this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index));

				index[0] = Ix; index[1] = Iy; index[2] = Iz + 1;
				m_RayIntersectionVoxels[2] =
					(this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index));

				index[0] = Ix; index[1] = Iy + 1; index[2] = Iz + 1;
				m_RayIntersectionVoxels[3] =
					(this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index));
			}
			else
			{
				m_RayIntersectionVoxels[0] =
					m_RayIntersectionVoxels[1] =
					m_RayIntersectionVoxels[2] =
					m_RayIntersectionVoxels[3] = nullptr;
			}
			break;
		}

		case 2:
		{
			if ((Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInx)
				&& (Iy >= 0) && (Iy < m_NumberOfVoxelsIny)
				&& (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInz))
			{
				index[0] = Ix; index[1] = Iy; index[2] = Iz;
				m_RayIntersectionVoxels[0] = (this->m_Image->GetBufferPointer()
					+ this->m_Image->ComputeOffset(index));

				index[0] = Ix + 1; index[1] = Iy; index[2] = Iz;
				m_RayIntersectionVoxels[1] = (this->m_Image->GetBufferPointer()
					+ this->m_Image->ComputeOffset(index));

				index[0] = Ix; index[1] = Iy; index[2] = Iz + 1;
				m_RayIntersectionVoxels[2] = (this->m_Image->GetBufferPointer()
					+ this->m_Image->ComputeOffset(index));

				index[0] = Ix + 1; index[1] = Iy; index[2] = Iz + 1;
				m_RayIntersectionVoxels[3] = (this->m_Image->GetBufferPointer()
					+ this->m_Image->ComputeOffset(index));
			}
			else
			{
				m_RayIntersectionVoxels[0] =
					m_RayIntersectionVoxels[1] =
					m_RayIntersectionVoxels[2] =
					m_RayIntersectionVoxels[3] = nullptr;
			}
			break;
		}

		case 3:
		{
			if ((Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInx)
				&& (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsIny)
				&& (Iz >= 0) && (Iz < m_NumberOfVoxelsInz))
			{
				index[0] = Ix; index[1] = Iy; index[2] = Iz;
				m_RayIntersectionVoxels[0] = (this->m_Image->GetBufferPointer()
					+ this->m_Image->ComputeOffset(index));

				index[0] = Ix + 1; index[1] = Iy; index[2] = Iz;
				m_RayIntersectionVoxels[1] = (this->m_Image->GetBufferPointer()
					+ this->m_Image->ComputeOffset(index));

				index[0] = Ix; index[1] = Iy + 1; index[2] = Iz;
				m_RayIntersectionVoxels[2] = (this->m_Image->GetBufferPointer()
					+ this->m_Image->ComputeOffset(index));

				index[0] = Ix + 1; index[1] = Iy + 1; index[2] = Iz;
				m_RayIntersectionVoxels[3] = (this->m_Image->GetBufferPointer()
					+ this->m_Image->ComputeOffset(index));
			}
			else
			{
				m_RayIntersectionVoxels[0] =
					m_RayIntersectionVoxels[1] =
					m_RayIntersectionVoxels[2] =
					m_RayIntersectionVoxels[3] = nullptr;
			}
			break;
		}
		}
	}
	//integrateabovethreshold

	double intensity;

	//  double posn3D_x, posn3D_y, posn3D_z;
	double sum;

	// Check if this is a valid ray
	if (!m_ValidRay)
	{
		return;
	}
	/* Step along the ray as quickly as possible
	   integrating the interpolated intensities. */

	for (m_NumVoxelPlanesTraversed = 0;
		m_NumVoxelPlanesTraversed < m_TotalRayVoxelPlanes;
		m_NumVoxelPlanesTraversed++)
	{
		double a, b, c, d;
		double y, z;

		a = (double)(*m_RayIntersectionVoxels[0]);
		b = (double)(*m_RayIntersectionVoxels[1] - a);
		c = (double)(*m_RayIntersectionVoxels[2] - a);
		d = (double)(*m_RayIntersectionVoxels[3] - a - b - c);

		switch (m_TraversalDirection)
		{
		case 1:
		{
			y = m_Position3Dvox[1].GetSum() - std::floor(m_Position3Dvox[1].GetSum());
			z = m_Position3Dvox[2].GetSum() - std::floor(m_Position3Dvox[2].GetSum());
			break;
		}
		case 2:
		{
			y = m_Position3Dvox[0].GetSum() - std::floor(m_Position3Dvox[0].GetSum());
			z = m_Position3Dvox[2].GetSum() - std::floor(m_Position3Dvox[2].GetSum());
			break;
		}
		case 3:
		{
			y = m_Position3Dvox[0].GetSum() - std::floor(m_Position3Dvox[0].GetSum());
			z = m_Position3Dvox[1].GetSum() - std::floor(m_Position3Dvox[1].GetSum());
			break;
		}
		}
		intensity = a + b * y + c * z + d * y*z;

		if (intensity > threshold)
		{
			sum += intensity - threshold;
		}
		double xBefore = m_Position3Dvox[0].GetSum();
		double yBefore = m_Position3Dvox[1].GetSum();
		double zBefore = m_Position3Dvox[2].GetSum();

		m_Position3Dvox[0] += m_VoxelIncrement[0];
		m_Position3Dvox[1] += m_VoxelIncrement[1];
		m_Position3Dvox[2] += m_VoxelIncrement[2];

		int dx = ((int)m_Position3Dvox[0].GetSum()) - ((int)xBefore);
		int dy = ((int)m_Position3Dvox[1].GetSum()) - ((int)yBefore);
		int dz = ((int)m_Position3Dvox[2].GetSum()) - ((int)zBefore);

		m_RayIntersectionVoxelIndex[0] += dx;
		m_RayIntersectionVoxelIndex[1] += dy;
		m_RayIntersectionVoxelIndex[2] += dz;

		int totalRayVoxelPlanes =
			dx + dy * m_NumberOfVoxelsInx + dz * m_NumberOfVoxelsInx * m_NumberOfVoxelsIny;

		m_RayIntersectionVoxels[0] += totalRayVoxelPlanes;
		m_RayIntersectionVoxels[1] += totalRayVoxelPlanes;
		m_RayIntersectionVoxels[2] += totalRayVoxelPlanes;
		m_RayIntersectionVoxels[3] += totalRayVoxelPlanes;
	}

	/* The ray passes through the volume one plane of voxels at a time,
	   however, if its moving diagonally the ray points will be further
	   apart so account for this by scaling by the distance moved. */

	pixval = sum.GetSum();
	pixval *= std::sqrt(m_VoxelIncrement[0] * spacing[0] * m_VoxelIncrement[0] * spacing[0]
		+ m_VoxelIncrement[1] * spacing[1] * m_VoxelIncrement[1] * spacing[1]
		+ m_VoxelIncrement[2] * spacing[2] * m_VoxelIncrement[2] * spacing[2]);


	//Assign the calculated value for the pixel to its corresponding position in the output array
	object2D[idx] = pixval;

}

float* calculateDRRwithCUDA(CUDAParamerters CUDA_Parameters, DRRParameters DRR_Parameters)
{

	cudaMemcpyToSymbol(d_DRR_Parameters, DRR_Parameters.stepInX, 12 * sizeof(float), 0, cudaMemcpyHostToDevice);

	//Block 6
	int num_Threads = CUDA_Parameters.numThreads;
	int num_Blocks = CUDA_Parameters.numBlocks;

	//------------------------------------------------------------
	//Launching the threads
	drrCUDA << < num_Blocks, num_Threads >> > (d_object2D, tex_object3D);
	//------------------------------------------------------------

	//Copying the result from the calculations from device to host
	long int vectorSize = (int)DRR_Parameters.size[0] * (int)DRR_Parameters.size[1];
	float* h_object2D = (float*)malloc(sizeof(float)*vectorSize);
	cudaMemcpy(h_object2D, d_object2D, vectorSize * sizeof(float), cudaMemcpyDeviceToHost);

	return h_object2D;
}


void HandleCudaKernelError(const cudaError_t CudaError, const char* pName /*= ""*/)
{
	if (CudaError == cudaSuccess)
		return;

	std::cerr << "The '" << pName << " kernel caused the following CUDA runtime error: " << cudaGetErrorString(CudaError) << std::endl;
}
