#include "drr.h"

drr::drr(std::string input)
{
	reader = nullptr;

	//定义像素类型，图像类型，三维有符号数，定义指针
	typedef float PixelType;
	const unsigned int Dimension = 3;
	typedef itk::Image< PixelType, Dimension > ImageType;
	typedef itk::ImageSeriesReader< ImageType > ReaderType;

	//声明读、写 DICOM 图 像 的 itk::GDCMImageIO对象
	//itk::GDCMSeriesFileNames对象将生成并将构成所有体数据的切片的文件名进行排序
	typedef itk::GDCMImageIO ImageIOType;
	typedef itk::GDCMSeriesFileNames NamesGeneratorType;
	ImageIOType::Pointer gdcmIO = ImageIOType::New();
	NamesGeneratorType::Pointer namesGenerator = NamesGeneratorType::New();

	bool result = reader.IsNull();
	double time = (double)clock() / CLOCKS_PER_SEC;
	//设置读取路径
	//用文件名发生器生成被读的文件名和被写的文件名
	namesGenerator->SetInputDirectory(input);
	const ReaderType::FileNamesContainer& filenames = namesGenerator->GetInputFileNames();

	//设置DICOM图像IO对象和被读的文件名的列表
	ReaderType::Pointer reader = ReaderType::New();
	reader->SetImageIO(gdcmIO);
	reader->SetFileNames(filenames);
	try
	{
		reader->Update();
	}
	catch (const itk::ExceptionObject & err)
	{
		std::string strerr = err.GetDescription();
		return;
	}

	image = reader->GetOutput();

	ImageType::SizeType sizeCT = image->GetLargestPossibleRegion().GetSize();
	InputImageType::SpacingType ctPixelSpacing = image->GetSpacing();
	//free(cpp_object3D);
	cpp_object3D = (float*)malloc(sizeof(float)*sizeCT[0] * sizeCT[1] * sizeCT[2]);
	for (int i = 0; i < sizeCT[0]; i++)
		for (int j = 0; j < sizeCT[1]; j++)
			for (int k = 0; k < sizeCT[2]; k++) {
				itk::Index<3> index;
				index[0] = i;
				index[1] = j;
				index[2] = k;
				cpp_object3D[k + j * (sizeCT[2]) + i * (sizeCT[1] * sizeCT[2])] = image->GetPixel(index);
			}

	int SizeCT[3];
	SizeCT[0] = sizeCT[0];
	SizeCT[1] = sizeCT[1];
	SizeCT[2] = sizeCT[2];
	float PixelSpacingCT[3];
	PixelSpacingCT[0] = ctPixelSpacing[0];
	PixelSpacingCT[1] = ctPixelSpacing[1];
	PixelSpacingCT[2] = ctPixelSpacing[2];
	loadDICOMInGPUMemory(cpp_object3D, SizeCT, PixelSpacingCT);
	free(cpp_object3D);
	cpp_object3D = NULL;
	std::cout << "readfile time:" << (double)clock() / CLOCKS_PER_SEC - time << "s" << std::endl;
}


drr::~drr()
{
	freeDICOMFromGPUMemory();
}

float* drr::cudaFloatDRR(float rx, float ry, float rz, int dx, int dy, float threshold)
{
	InputImageType::PointType ori;
	ori[0] = 0.0;
	ori[1] = 0.0;
	ori[2] = 0.0;
	image->SetOrigin(ori);

	typedef itk::Euler3DTransform<double> TransformType;
	TransformType::Pointer transform = TransformType::New();
	TransformType::OutputVectorType translation;
	translation[0] = 0;
	translation[1] = 0;
	translation[2] = 0;
	const double dtr = (atan(1.0)*4.0) / 180.0;
	transform->SetTranslation(translation);
	//transform->SetComputeZYX(true);
	transform->SetRotation(dtr*rx, dtr*ry, dtr*rz);

	itk::Vector<double> SourceWorld;

	typename InputImageType::SizeType sizeCT;
	typename InputImageType::RegionType regionCT;
	typename InputImageType::SpacingType ctPixelSpacing;
	typename InputImageType::PointType ctOrigin;

	ctPixelSpacing = image->GetSpacing();
	ctOrigin = image->GetOrigin();
	regionCT = image->GetLargestPossibleRegion();
	sizeCT = regionCT.GetSize();
	InputImageType::SizeType imSize = image->GetBufferedRegion().GetSize();

	TransformType::InputPointType isocenter;
	isocenter[0] = ctOrigin[0] + ctPixelSpacing[0] * static_cast<double>(imSize[0]) / 2.0;
	isocenter[1] = ctOrigin[1] + ctPixelSpacing[1] * static_cast<double>(imSize[1]) / 2.0;
	isocenter[2] = ctOrigin[2] + ctPixelSpacing[2] * static_cast<double>(imSize[2]) / 2.0;
	transform->SetCenter(isocenter);

	TransformType::Pointer m_InverseTransform = TransformType::New();
	m_InverseTransform->SetComputeZYX(true);
	TransformType::Pointer m_ComposeTransform = TransformType::New();
	m_ComposeTransform->SetComputeZYX(true);
	TransformType::Pointer m_GantryRoTransform = TransformType::New();
	m_GantryRoTransform->SetComputeZYX(true);
	m_GantryRoTransform->SetIdentity();
	TransformType::Pointer m_CamShiftTransform = TransformType::New();
	m_CamShiftTransform->SetComputeZYX(true);
	m_CamShiftTransform->SetIdentity();
	TransformType::Pointer m_CamRotTransform = TransformType::New();
	m_CamRotTransform->SetComputeZYX(true);
	m_CamRotTransform->SetIdentity();
	float dtrf = (atan(1.0)*4.0) / 180.0;
	m_CamRotTransform->SetRotation(dtrf*(-90.0), 0.0, 0.0);

	m_ComposeTransform->SetIdentity();
	m_ComposeTransform->Compose(transform, 0);
	float rprojection = 90.;
	float angle = -dtr * rprojection;
	m_GantryRoTransform->SetRotation(angle, -angle, 0.);
	m_GantryRoTransform->SetCenter(isocenter);
	m_ComposeTransform->Compose(m_GantryRoTransform, 0);
	TransformType::OutputVectorType focalpointtranslation;
	focalpointtranslation[0] = -isocenter[0];
	focalpointtranslation[1] = 400.0 - isocenter[1];
	focalpointtranslation[2] = -isocenter[2];
	//std::cout << isocenter[0] << " " << isocenter[1] << " " << isocenter[2] << endl;
	//std::cout << focalpointtranslation[0] << " " << focalpointtranslation[1] << " " << focalpointtranslation[2] << endl;
	m_CamShiftTransform->SetTranslation(focalpointtranslation);
	m_ComposeTransform->Compose(m_CamShiftTransform, 0);
	m_ComposeTransform->Compose(m_CamRotTransform, 0);



	m_ComposeTransform->GetInverse(m_InverseTransform);

	itk::Matrix<double> m_matrix = m_InverseTransform->GetMatrix();
	itk::Vector<double,3> m_offset= m_InverseTransform->GetOffset();
	
	PointType m_SourcePoint;
	m_SourcePoint[0] = 0.;
	m_SourcePoint[1] = 0.;
	m_SourcePoint[2] = 0.;
	PointType tmp_SourceWorld;
	m_InverseTransform->TransformPoint(m_SourcePoint);
	tmp_SourceWorld = m_matrix * m_SourcePoint + m_offset;

	itk::Matrix<double> m_direction = image->GetDirection();

	//output_2d_image
	double output_spacing[3];
	output_spacing[0] = 0.3125;//sx
	output_spacing[1] = 0.3125;//sy
	output_spacing[2] = 1;
	double output_size[3];
	output_size[0] = dx;//dx
	output_size[1] = dy;//dy
	output_size[2] = 1.0;

	itk::Vector<double> origin;
	origin[0] = -output_spacing[0] * (output_size[0] - 1.) / 2.;
	origin[1] = -output_spacing[1] * (output_size[1] - 1.) / 2.;
	origin[2] = -400;

	itk::Vector<double> corner00;
	itk::Vector<double> cornerOrigin01;
	itk::Vector<double> corner01;
	itk::Vector<double> cornerOrigin10;
	itk::Vector<double> corner10;

	cornerOrigin01[0] = origin[0] + output_spacing[0] * (output_size[0] - 1.);
	cornerOrigin01[1] = origin[1];
	cornerOrigin01[2] = origin[2];

	cornerOrigin10[0] = origin[0];
	cornerOrigin10[1] = origin[1] + output_spacing[1] * (output_size[1] - 1.);
	cornerOrigin10[2] = origin[2];

	corner00 = m_matrix * origin + m_offset;
	corner01 = m_matrix * cornerOrigin01 + m_offset;
	corner10 = m_matrix * cornerOrigin10 + m_offset;
	SourceWorld[0] = tmp_SourceWorld[0];
	SourceWorld[1] = tmp_SourceWorld[1];
	SourceWorld[2] = tmp_SourceWorld[2];

	itk::Vector<double> stepInX;
	stepInX[0] = (corner01[0] - corner00[0]) / (output_size[0]);
	stepInX[1] = (corner01[1] - corner00[1]) / (output_size[0]);
	stepInX[2] = (corner01[2] - corner00[2]) / (output_size[0]);
	itk::Vector<double> stepInY;
	stepInY[0] = (corner10[0] - corner00[0]) / (output_size[1]);
	stepInY[1] = (corner10[1] - corner00[1]) / (output_size[1]);
	stepInY[2] = (corner10[2] - corner00[2]) / (output_size[1]);

	DRRParameters para;
	para.stepInX[0] = stepInX[0];
	para.stepInX[1] = stepInX[1];
	para.stepInX[2] = stepInX[2];
	para.stepInY[0] = stepInY[0];
	para.stepInY[1] = stepInY[1];
	para.stepInY[2] = stepInY[2];
	para.corner00[0] = corner00[0];
	para.corner00[1] = corner00[1];
	para.corner00[2] = corner00[2];
	para.SourceWorld[0] = SourceWorld[0];
	para.SourceWorld[1] = SourceWorld[1];
	para.SourceWorld[2] = SourceWorld[2];
	para.size[0] = output_size[0];
	para.size[1] = output_size[1];
	para.threshold = threshold;

	
	Image3D image3d;
	image3d.image = cpp_object3D;
	image3d.isoCenter[0] = isocenter[0];
	image3d.isoCenter[1] = isocenter[1];
	image3d.isoCenter[2] = isocenter[2];
	image3d.PixelSpacingCT[0] = ctPixelSpacing[0];
	image3d.PixelSpacingCT[1] = ctPixelSpacing[1];
	image3d.PixelSpacingCT[2] = ctPixelSpacing[2];
	image3d.SizeCT[0] = sizeCT[0];
	image3d.SizeCT[1] = sizeCT[1];
	image3d.SizeCT[2] = sizeCT[2];

	CUDAParamerters cudaPara;
	cudaDeviceProp prop;
	cudaGetDeviceProperties(&prop, 0);
	cudaPara.numThreads = prop.maxThreadsPerBlock;
	cudaPara.numBlocks = (int)ceil((float)output_size[0] * output_size[1] / cudaPara.numThreads);


	loadOuputVariablesInGPUMemory((int)output_size[0],(int)output_size[1]);
	float* object2d = calFloatDRRwithCUDA(cudaPara, para);


	freeAuxiliaryVariablesInGPUMemory();

	return object2d;

	//////////////以下代码为保存DRR到硬盘
	/*
	typedef itk::Image< unsigned char, 2 > ImageType;
	ImageType::Pointer image2d = ImageType::New();

	ImageType::IndexType start;
	start[0] = 0;  // first index on X
	start[1] = 0;  // first index on Y

	ImageType::SizeType  size;
	size[0] = 1024;  // size along X
	size[1] = 1024;  // size along Y

	ImageType::RegionType region;
	region.SetSize(size);
	region.SetIndex(start);
	image2d->SetRegions(region);
	image2d->Allocate();

	ImageType::PixelType  initialValue = 0;
	image2d->FillBuffer(initialValue);

	int k = 0;
	for (int i = 0; i < size[1]; i++)
	{
		for (int j = size[0]-1; j >= 0; j--)
		{
			ImageType::IndexType pixelIndex;
			pixelIndex[0] = i;
			pixelIndex[1] = j;
			image2d->SetPixel(pixelIndex, object2d[k]);
			k++;
		}
	}
	

	itk::BMPImageIOFactory::RegisterOneFactory();
	typedef itk::ImageFileWriter< ImageType > WriterType;
	WriterType::Pointer writer = WriterType::New();

	// Or other file format, such as mhd/mha, png, jpg etc.
	writer->SetFileName("output.bmp");
	writer->SetInput(image2d);
	writer->Update();
	*/
	

}

unsigned char * drr::cudaUCharDRR(float rx, float ry, float rz, int dx, int dy, float threshold)
{
	InputImageType::PointType ori;
	ori[0] = 0.0;
	ori[1] = 0.0;
	ori[2] = 0.0;
	image->SetOrigin(ori);

	typedef itk::Euler3DTransform<double> TransformType;
	TransformType::Pointer transform = TransformType::New();
	TransformType::OutputVectorType translation;
	translation[0] = 0;
	translation[1] = 0;
	translation[2] = 0;
	const double dtr = (atan(1.0)*4.0) / 180.0;
	transform->SetTranslation(translation);
	//transform->SetComputeZYX(true);
	transform->SetRotation(dtr*rx, dtr*ry, dtr*rz);

	itk::Vector<double> SourceWorld;

	typename InputImageType::SizeType sizeCT;
	typename InputImageType::RegionType regionCT;
	typename InputImageType::SpacingType ctPixelSpacing;
	typename InputImageType::PointType ctOrigin;

	ctPixelSpacing = image->GetSpacing();
	ctOrigin = image->GetOrigin();
	regionCT = image->GetLargestPossibleRegion();
	sizeCT = regionCT.GetSize();
	InputImageType::SizeType imSize = image->GetBufferedRegion().GetSize();

	TransformType::InputPointType isocenter;
	isocenter[0] = ctOrigin[0] + ctPixelSpacing[0] * static_cast<double>(imSize[0]) / 2.0;
	isocenter[1] = ctOrigin[1] + ctPixelSpacing[1] * static_cast<double>(imSize[1]) / 2.0;
	isocenter[2] = ctOrigin[2] + ctPixelSpacing[2] * static_cast<double>(imSize[2]) / 2.0;
	transform->SetCenter(isocenter);

	TransformType::Pointer m_InverseTransform = TransformType::New();
	m_InverseTransform->SetComputeZYX(true);
	TransformType::Pointer m_ComposeTransform = TransformType::New();
	m_ComposeTransform->SetComputeZYX(true);
	TransformType::Pointer m_GantryRoTransform = TransformType::New();
	m_GantryRoTransform->SetComputeZYX(true);
	m_GantryRoTransform->SetIdentity();
	TransformType::Pointer m_CamShiftTransform = TransformType::New();
	m_CamShiftTransform->SetComputeZYX(true);
	m_CamShiftTransform->SetIdentity();
	TransformType::Pointer m_CamRotTransform = TransformType::New();
	m_CamRotTransform->SetComputeZYX(true);
	m_CamRotTransform->SetIdentity();
	float dtrf = (atan(1.0)*4.0) / 180.0;
	m_CamRotTransform->SetRotation(dtrf*(-90.0), 0.0, 0.0);

	m_ComposeTransform->SetIdentity();
	m_ComposeTransform->Compose(transform, 0);
	float rprojection = 90.;
	float angle = -dtr * rprojection;
	m_GantryRoTransform->SetRotation(angle, 2*angle, 0.);
	m_GantryRoTransform->SetCenter(isocenter);
	m_ComposeTransform->Compose(m_GantryRoTransform, 0);
	TransformType::OutputVectorType focalpointtranslation;
	focalpointtranslation[0] = -isocenter[0];
	focalpointtranslation[1] = 400.0 - isocenter[1];
	focalpointtranslation[2] = -isocenter[2];
	//std::cout << isocenter[0] << " " << isocenter[1] << " " << isocenter[2] << endl;
	//std::cout << focalpointtranslation[0] << " " << focalpointtranslation[1] << " " << focalpointtranslation[2] << endl;
	m_CamShiftTransform->SetTranslation(focalpointtranslation);
	m_ComposeTransform->Compose(m_CamShiftTransform, 0);
	m_ComposeTransform->Compose(m_CamRotTransform, 0);



	m_ComposeTransform->GetInverse(m_InverseTransform);

	itk::Matrix<double> m_matrix = m_InverseTransform->GetMatrix();
	itk::Vector<double, 3> m_offset = m_InverseTransform->GetOffset();

	PointType m_SourcePoint;
	m_SourcePoint[0] = 0.;
	m_SourcePoint[1] = 0.;
	m_SourcePoint[2] = 0.;
	PointType tmp_SourceWorld;
	m_InverseTransform->TransformPoint(m_SourcePoint);
	tmp_SourceWorld = m_matrix * m_SourcePoint + m_offset;

	itk::Matrix<double> m_direction = image->GetDirection();

	//output_2d_image
	double output_spacing[3];
	output_spacing[0] = 0.3125;//sx
	output_spacing[1] = 0.3125;//sy
	output_spacing[2] = 1;
	double output_size[3];
	output_size[0] = dx;//dx
	output_size[1] = dy;//dy
	output_size[2] = 1.0;

	itk::Vector<double> origin;
	origin[0] = -output_spacing[0] * (output_size[0] - 1.) / 2.;
	origin[1] = -output_spacing[1] * (output_size[1] - 1.) / 2.;
	origin[2] = -400;

	itk::Vector<double> corner00;
	itk::Vector<double> cornerOrigin01;
	itk::Vector<double> corner01;
	itk::Vector<double> cornerOrigin10;
	itk::Vector<double> corner10;

	cornerOrigin01[0] = origin[0] + output_spacing[0] * (output_size[0] - 1.);
	cornerOrigin01[1] = origin[1];
	cornerOrigin01[2] = origin[2];

	cornerOrigin10[0] = origin[0];
	cornerOrigin10[1] = origin[1] + output_spacing[1] * (output_size[1] - 1.);
	cornerOrigin10[2] = origin[2];

	corner00 = m_matrix * origin + m_offset;
	corner01 = m_matrix * cornerOrigin01 + m_offset;
	corner10 = m_matrix * cornerOrigin10 + m_offset;
	SourceWorld[0] = tmp_SourceWorld[0];
	SourceWorld[1] = tmp_SourceWorld[1];
	SourceWorld[2] = tmp_SourceWorld[2];

	itk::Vector<double> stepInX;
	stepInX[0] = (corner01[0] - corner00[0]) / (output_size[0]);
	stepInX[1] = (corner01[1] - corner00[1]) / (output_size[0]);
	stepInX[2] = (corner01[2] - corner00[2]) / (output_size[0]);
	itk::Vector<double> stepInY;
	stepInY[0] = (corner10[0] - corner00[0]) / (output_size[1]);
	stepInY[1] = (corner10[1] - corner00[1]) / (output_size[1]);
	stepInY[2] = (corner10[2] - corner00[2]) / (output_size[1]);

	DRRParameters para;
	para.stepInX[0] = stepInX[0];
	para.stepInX[1] = stepInX[1];
	para.stepInX[2] = stepInX[2];
	para.stepInY[0] = stepInY[0];
	para.stepInY[1] = stepInY[1];
	para.stepInY[2] = stepInY[2];
	para.corner00[0] = corner00[0];
	para.corner00[1] = corner00[1];
	para.corner00[2] = corner00[2];
	para.SourceWorld[0] = SourceWorld[0];
	para.SourceWorld[1] = SourceWorld[1];
	para.SourceWorld[2] = SourceWorld[2];
	para.size[0] = output_size[0];
	para.size[1] = output_size[1];
	para.threshold = threshold;


	Image3D image3d;
	image3d.image = cpp_object3D;
	image3d.isoCenter[0] = isocenter[0];
	image3d.isoCenter[1] = isocenter[1];
	image3d.isoCenter[2] = isocenter[2];
	image3d.PixelSpacingCT[0] = ctPixelSpacing[0];
	image3d.PixelSpacingCT[1] = ctPixelSpacing[1];
	image3d.PixelSpacingCT[2] = ctPixelSpacing[2];
	image3d.SizeCT[0] = sizeCT[0];
	image3d.SizeCT[1] = sizeCT[1];
	image3d.SizeCT[2] = sizeCT[2];

	CUDAParamerters cudaPara;
	cudaDeviceProp prop;
	cudaGetDeviceProperties(&prop, 0);
	cudaPara.numThreads = prop.maxThreadsPerBlock;
	cudaPara.numBlocks = (int)ceil((float)output_size[0] * output_size[1] / cudaPara.numThreads);


	loadOuputVariablesInGPUMemory((int)output_size[0], (int)output_size[1], 1);
	unsigned char* object2d = calUCharDRRwithCUDA(cudaPara, para);


	freeAuxiliaryVariablesInGPUMemory(1);

	return object2d;
}
