ctconv_reconstuct.cpp

一个以研究CT重建算法为主的仿真演示程序 C++命令行形式 详见readme thanks

			{
				double maxproj=projfunc[0];
				double minproj=projfunc[0];

				for (i=0;i<width;i++)
				{
					if (maxproj<projfunc[i]) maxproj=projfunc[i];
					if (minproj>projfunc[i]) minproj=projfunc[i];
				}

				for (i=0;i<width;i++)
				{
					projfunc[i]=projfunc[i]+rand()/(double)RAND_MAX*maxproj/SNR;
				}
			}
			//apply the R-L convolution
			//the R-L convolution is symmetric,no need to converse
			double projtemp=0;
			double projfuncRL[width+2*RLlen];	//create a projection function array to store the values after RL convolution

			for(i=0;i<width+2*RLlen;i++)
			{
				if(i<RLlen) projfuncRL[i]=0.;
				else if(i>=(width+RLlen)) projfuncRL[i]=0.;
				else projfuncRL[i]=projfunc[i-RLlen];
			}

			for(i=RLlen;i<RLlen+width;i++)
			{
				for(j=0;j<2*RLlen+1;j++)
				{
					projtemp=projtemp+projfuncRL[i+j-RLlen]*rlconv[j];
				}
				projfunc[i-RLlen]=projtemp;
				projtemp=0;
			}

		}


		if (ctrl=='2')				//S-L convolution no noises
		{
			//apply the S-L convolution
			if (ctrln=='y')	//there is noises
			{
				double maxproj=projfunc[0];
				double minproj=projfunc[0];

				for (i=0;i<width;i++)
				{
					if (maxproj<projfunc[i]) maxproj=projfunc[i];
					if (minproj>projfunc[i]) minproj=projfunc[i];
				}

				for (i=0;i<width;i++)
				{
					projfunc[i]=projfunc[i]+rand()/(double)RAND_MAX*maxproj/SNR;
				}
			}

			double projtemp=0.;
			double projfuncSL[width+2*SLlen];	//create a projection function array to store the values after SL convolution

			for(i=0;i<width+2*SLlen;i++)
			{
				if(i<SLlen) projfuncSL[i]=0.;
				else if(i>=(width+SLlen)) projfuncSL[i]=0.;
				else projfuncSL[i]=projfunc[i-SLlen];
			}

			for(i=SLlen;i<SLlen+width;i++)
			{
				for(j=0;j<2*SLlen+1;j++)
				{
					projtemp=projtemp+projfuncSL[i+j-SLlen]*slconv[j];
				}
				projfunc[i-SLlen]=projtemp;
				projtemp=0.;
			}
		}

		if (ctrl=='3')				//S-L convolution no noises
		{
			//apply the S-L convolution
			if (ctrln=='y')	//there is noises
			{
				double maxproj=projfunc[0];
				double minproj=projfunc[0];

				for (i=0;i<width;i++)
				{
					if (maxproj<projfunc[i]) maxproj=projfunc[i];
					if (minproj>projfunc[i]) minproj=projfunc[i];
				}

				for (i=0;i<width;i++)
				{
					projfunc[i]=projfunc[i]+rand()/(double)RAND_MAX*maxproj/SNR;
				}
			}

			double projtemp=0.;
			double projfuncmyconv[width+2*MYCONVLEN];	//create a projection function array to store the values after SL convolution

			for(i=0;i<width+2*MYCONVLEN;i++)
			{
				if(i<MYCONVLEN) projfuncmyconv[i]=0.;
				else if(i>=(width+MYCONVLEN)) projfuncmyconv[i]=0.;
				else projfuncmyconv[i]=projfunc[i-MYCONVLEN];
			}

			for(i=MYCONVLEN;i<MYCONVLEN+width;i++)
			{
				if (ctrlmy=='0')
					for(j=0;j<2*MYCONVLEN+1;j++) projtemp=projtemp+projfuncmyconv[i+j-MYCONVLEN]*convsharphighpass[j];
				if (ctrlmy=='1')
					for(j=0;j<2*MYCONVLEN+1;j++) projtemp=projtemp+projfuncmyconv[i+j-MYCONVLEN]*convsmoothhighpass[j];
				if (ctrlmy=='2')
					for(j=0;j<2*MYCONVLEN+1;j++) projtemp=projtemp+projfuncmyconv[i+j-MYCONVLEN]*convalteredhighpass[j];

				projfunc[i-MYCONVLEN]=projtemp;
				projtemp=0.;
			}
		}

		if (ctrl=='0')
		{
			if (ctrln=='y')	//there is noises
			{
				double maxproj=projfunc[0];
				double minproj=projfunc[0];

				for (i=0;i<width;i++)
				{
					if (maxproj<projfunc[i]) maxproj=projfunc[i];
					if (minproj>projfunc[i]) minproj=projfunc[i];
				}

				for (i=0;i<width;i++)
				{
					projfunc[i]=projfunc[i]+rand()/(double)RAND_MAX*maxproj/SNR;
				}
			}

		}

		double maxproj=projfunc[0];
		double minproj=projfunc[0];

		for (i=0;i<width;i++)
		{
			if (maxproj<projfunc[i]) maxproj=projfunc[i];
			if (minproj>projfunc[i]) minproj=projfunc[i];
		}

		double diffproj=maxproj-minproj;
		for (i=0;i<width;i++)
		{
			projfunc[i]=(projfunc[i]-minproj)/diffproj*255;
		}		

		//back projection to temp image
		for (j=0;j<wt;j++)	
			for(i=0;i<ht;i++)
				for(k=0;k<cht;k++)
					datat[i*st+j*cht+k]=projfunc[j];		//one rotated image was projected! stored in temp

		//now rotate back
		float mb[6];
		CvMat Mb = cvMat (2, 3, CV_32F, mb);	  
		if (opt)		// 旋转加缩放	    
			factor = (cos (angle * CV_PI / 180.) + 1.0) * 2;	  
		else			//  仅仅旋转	     
			factor = 1;	  
		mb[0] = (float) (factor * cos (-angle * CV_PI / 180.));	//clockwise -angle degrees	  
		mb[1] = (float) (factor * sin (-angle * CV_PI / 180.));	  
		mb[3] = -mb[1];	  
		mb[4] = mb[0];	  
		// 将旋转中心移至图像中间	  
		mb[2] = wr * 0.5f;	  
		mb[5] = hr * 0.5f;	  //  dst(x,y) = A * src(x,y) + b
		cvZero (temp2);
		cvGetQuadrangleSubPix (temp,temp2,&Mb);
		
		//calculate the projection image
		for (i=0;i<hp1;i++)	
			for(j=0;j<wp1;j++)
				for(k=0;k<chp1;k++)
					datap1temp[i*sp1+j*chp1+k]=datap1temp[i*sp1+j*chp1+k]+datat2[i*st+j*cht+k];
		
		//the angle increases
		angle = (int) (angle + delta) % 360;
	}
	

	/*for (i=0;i<hp1;i++)	// all over and process the final image
		for(j=0;j<wp1;j++)
			for(k=0;k<chp1;k++)
				datap1[i*sp1+j*chp1+k]=datap1[i*sp1+j*chp1+k]/times;*/

	//how to handle this problem:the processing of the pre-finished image
	double min=datap1temp[0];
	double max=datap1temp[0];
	for (i=0;i<hp1;i++)	// all over and process the final image
	{
		for(j=0;j<wp1;j++)
		{
			for(k=0;k<chp1;k++)
			{
				if (min>datap1temp[i*sp1+j*chp1+k])
					min=datap1temp[i*sp1+j*chp1+k];
				if (max<datap1temp[i*sp1+j*chp1+k])
					max=datap1temp[i*sp1+j*chp1+k];
			}
		}
	}

	double diff=max-min;
	for (i=0;i<hp1;i++)	// extend the grey value
		for(j=0;j<wp1;j++)
			for(k=0;k<chp1;k++)
				datap1[i*sp1+j*chp1+k]=(datap1temp[i*sp1+j*chp1+k]-min)/diff*255;

	//text="My Conv Back Projection Reconstructs S-L Head Model";
	//cvPutText(img1proj1, text, cvPoint(10,20), &myfont, fontcolor );

	finish=clock();
	totaltime=(double)(finish-start)/CLOCKS_PER_SEC;
	cout<<"\nThe programme consumes "<<totaltime<<"s"<<endl;

	cout<<"Press any key to destroy the windows and realease the memory"<<endl;

	cvNamedWindow("Projection function image",1);
	cvMoveWindow("Projection function image",200,100);
	cvShowImage("Projection function image",temp2);


	//create gui window and display the original image
	cvNamedWindow("Original Image of S-L Head model",1);
	cvMoveWindow("Original Image of S-L Head model",100,100);
	cvShowImage("Original Image of S-L Head model",img1);
	//save the original image
	cvSaveImage(saveimg1,img1);


	cvNamedWindow("Back Projection Image",1);
	cvMoveWindow("Back Projection Image",300,100);
	cvShowImage("Back Projection Image",img1proj1);
	//save the projection image
	cvSaveImage(saveimg1proj1,img1proj1);
	
	cvWaitKey(0);

	cvDestroyWindow("Back Projection Image");
	cvReleaseImage(&img1proj1);

	cvDestroyWindow("Projection function image");
	cvReleaseImage(&temp2);

	cvDestroyWindow("Original Image of S-L Head model");
	cvReleaseImage(&img1);

	cvReleaseImage(&temp);

	delete(datap1temp);

	cout<<endl<<"Press any other key to continue and press 'q' to quit"<<endl;
	cin>>ctrlquit;
	}while (ctrlquit!='q');

	return 0;
}