ctconv_reconstuct.cpp

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

	//define the rotation angle of the ellipse objects
	double alphaa=90.,alphab=90.,alphac=72.,alphad=108.,alphae=90.,alphaf=0.,alphag=0.,alphah=0.,alphai=0.,alphaj=90.;
	//define the size of the image
		
	//the original image
	IplImage* img1;
	img1=cvCreateImage(cvSize(width,height),imagedepth,1);
	//get the information
	int h=img1->height;
	int w=img1->width;
	int s=img1->widthStep;
	int ch=img1->nChannels;
	unsigned char *data;
	data=(unsigned char *)img1->imageData;

	//initialize the backround if Background is larger,the picture is lighter
	for (i=0;i<h;i++)
		for(j=0;j<w;j++)
			for(k=0;k<ch;k++)
				data[i*s+j*ch+k]=Background;

	//plot some ellipse according to the model
	cvEllipse(img1,pta,sizea,alphaa,0,360,colora,-1,imagedepth,0);
	cvEllipse(img1,ptb,sizeb,alphab,0,360,colorb,-1,imagedepth,0);
	cvEllipse(img1,ptc,sizec,alphac,0,360,colorc,-1,imagedepth,0);
	cvEllipse(img1,ptd,sized,alphad,0,360,colord,-1,imagedepth,0);
	cvEllipse(img1,pte,sizee,alphae,0,360,colore,-1,imagedepth,0);
	cvEllipse(img1,ptf,sizef,alphaf,0,360,colorf,-1,imagedepth,0);
	cvEllipse(img1,ptg,sizeg,alphag,0,360,colorg,-1,imagedepth,0);
	cvEllipse(img1,pth,sizeh,alphah,0,360,colorh,-1,imagedepth,0);
	cvEllipse(img1,pti,sizei,alphai,0,360,colori,-1,imagedepth,0);
	cvEllipse(img1,ptj,sizej,alphaj,0,360,colorj,-1,imagedepth,0);

	double hscale=1.5;
	double vscale=1.5;
	double shear=0.2;
	int thickness=1;
	int line_type=8;

	CvFont myfont;
	cvInitFont(&myfont,CV_FONT_HERSHEY_PLAIN,hscale,vscale,shear,thickness,line_type);
	
	//put the comment text upon the image
	char* text="Original Image of S-L Head Model----Plotted by Simon";
	CvScalar fontcolor=cvScalar(int(pow(2.,imagedepth)-1));
	cvPutText(img1, text, cvPoint(50,20), &myfont, fontcolor );

	//now the plotting is over

	//now select the reconstruct methods
	cout<<endl<<"Now please select your reconstruction method:"<<endl<<"0-->Direct back projection"<<endl<<"1-->R-L convolution"
		<<endl<<"2-->S-L convolution"<<endl<<"3-->My high pass convolution"<<endl;
	do{
		cin>>ctrl;
		if ((ctrl!='0')&&(ctrl!='1')&&(ctrl!='2')&&(ctrl!='3')) cout<<endl<<"Invalid Command! Please Renter:";
	}while((ctrl!='0')&&(ctrl!='1')&&(ctrl!='2')&&(ctrl!='3'));

	if (ctrl=='3')
	{
		cout<<endl<<"Now Please select the type of the high pass convolution"<<endl<<"0-->Sharp high pass convolution"<<endl
			<<"1-->Smooth high pass convolution"<<endl<<"2-->Altered high pass convolution"<<endl;
		do{
		cin>>ctrlmy;
		if ((ctrlmy!='0')&&(ctrlmy!='1')&&(ctrlmy!='2')) cout<<endl<<"Invalid Command! Please Renter:";
		}while ((ctrlmy!='0')&&(ctrlmy!='1')&&(ctrlmy!='2'));
	}


	cout<<endl<<"Now please select whether noises or not? [y/n]:";
	do{
		cin>>ctrln;
		if ((ctrln!='y')&&(ctrln!='n')) cout<<endl<<"Invalid Command! Please Renter:";
	}while((ctrln!='y')&&(ctrln!='n'));

	if (ctrln=='y')
	{
		cout<<endl<<"Please enter the SNR of the noises(larger than 10 recommended):";
		do{
			cin>>SNR;
			if (SNR<0) cout<<"Invalid Command! Please Renter:";
		}while(SNR<0);
	}

	//create the direct back projection image
	IplImage* img1proj1;
	img1proj1=cvCreateImage(cvSize(width,height),imagedepth,1);

	int hp1=img1proj1->height;
	int wp1=img1proj1->width;
	int sp1=img1proj1->widthStep;
	int chp1=img1proj1->nChannels;
	unsigned char *datap1;
	datap1=(unsigned char *)img1proj1->imageData;

	//initialize the projection image
	for (i=0;i<hp1;i++)	
		for(j=0;j<wp1;j++)
			for(k=0;k<chp1;k++)
				datap1[i*sp1+j*chp1+k]=0;


	//create a array as the projection function
	double projfunc[height];
	for (i=0;i<hp1;i++)	
		projfunc[i]=0;		//initialize the projection function


	//rotate the image to get the different projection function
	//create the rotated image
	IplImage* img1rotate1 = 0;
	img1rotate1=cvCreateImage(cvSize(width,height),imagedepth,1);
	
	unsigned int delta;
	cout<<endl<<"Please enter the scan step angle:";
	do{
		cin>>delta;
		if ((delta<=0)||(delta>360)) cout<<endl<<"Invalid Command! Please Renter:";
	}while ((delta<=0)||(delta>360));

	int angle = 0;      
	int opt = 0;		// 1： 旋转加缩放      // 0:  仅仅旋转      
	double factor;       

	int wr = img1rotate1->width;	  
	int hr = img1rotate1->height;
	int sr = img1rotate1->widthStep;
	int chr = img1rotate1->nChannels;
	unsigned char *rotatedata;
	rotatedata=(unsigned char *)img1rotate1->imageData;

	//create a temp image to record the projection function on a image when each rotation takes place
	IplImage* temp = 0;
	temp=cvCreateImage(cvSize(width,height),imagedepth,1);
	
	int ht=temp->height;
	int wt=temp->width;
	int st=temp->widthStep;
	int cht=temp->nChannels;
	unsigned char *datat;
	datat=(unsigned char *)temp->imageData;
	
	//create a temp2 image to record temp when rotated back
	IplImage* temp2 = 0;
	temp2=cvCreateImage(cvSize(width,height),imagedepth,1);

	int ht2=temp2->height;
	int wt2=temp2->width;
	int st2=temp2->widthStep;
	int cht2=temp2->nChannels;
	unsigned char *datat2;
	datat2=(unsigned char *)temp2->imageData;

	//initialize the temp image
	for (i=0;i<ht;i++)	
		for(j=0;j<wt;j++)
			for(k=0;k<cht;k++)
				datat[i*st+j*cht+k]=0;

	//because the datap1 automatically change the data between,need a temp array to record for afterword final processing
	double *datap1temp;
	datap1temp=new double[N];
	for (i=0;i<N;i++)
		datap1temp[i]=0.;	//initialize to 0


	cout<<endl<<"Processing,please stand by...";

	clock_t start,finish;
	start=clock();
	
	int times=360/delta;	//rotate how many times
	for (t=0;t<times;t++)	
	{
		for (i=0;i<height;i++)	//initialize the projection function
			projfunc[i]=0;

		for (i=0;i<ht;i++)	//initialize the temp image
			for(j=0;j<wt;j++)
				for(k=0;k<cht;k++)
					datat[i*st+j*cht+k]=0;
		

		//calculate the projection function
		for (i=0;i<height;i++)
		{
			double ba=pow((sizea.height*sin((angle-alphaa)*CV_PI/180.)),2);
			double aa=pow((sizea.width*cos((angle-alphaa)*CV_PI/180.)),2);
			double ca=pow(((i-400)-(pta.x-400)*cos(angle*CV_PI/180.)-(400-pta.y)*sin(angle*CV_PI/180.)),2);

			double bb=pow((sizeb.height*sin((angle-alphab)*CV_PI/180.)),2);
			double ab=pow((sizeb.width*cos((angle-alphab)*CV_PI/180.)),2);
			double cb=pow(((i-400)-(ptb.x-400)*cos(angle*CV_PI/180.)-(400-ptb.y)*sin(angle*CV_PI/180.)),2);

			double bc=pow((sizec.height*sin((angle-alphac)*CV_PI/180.)),2);
			double ac=pow((sizec.width*cos((angle-alphac)*CV_PI/180.)),2);
			double cc=pow(((i-400)-(ptc.x-400)*cos(angle*CV_PI/180.)-(400-ptc.y)*sin(angle*CV_PI/180.)),2);

			double bd=pow((sized.height*sin((angle-alphad)*CV_PI/180.)),2);
			double ad=pow((sized.width*cos((angle-alphad)*CV_PI/180.)),2);
			double cd=pow(((i-400)-(ptd.x-400)*cos(angle*CV_PI/180.)-(400-ptd.y)*sin(angle*CV_PI/180.)),2);

			double be=pow((sizee.height*sin((angle-alphae)*CV_PI/180.)),2);
			double ae=pow((sizee.width*cos((angle-alphae)*CV_PI/180.)),2);
			double ce=pow(((i-400)-(pte.x-400)*cos(angle*CV_PI/180.)-(400-pte.y)*sin(angle*CV_PI/180.)),2);

			double bf=pow((sizef.height*sin((angle-alphaf)*CV_PI/180.)),2);
			double af=pow((sizef.width*cos((angle-alphaf)*CV_PI/180.)),2);
			double cf=pow(((i-400)-(ptf.x-400)*cos(angle*CV_PI/180.)-(400-ptf.y)*sin(angle*CV_PI/180.)),2);

			double bg=pow((sizeg.height*sin((angle-alphag)*CV_PI/180.)),2);
			double ag=pow((sizeg.width*cos((angle-alphag)*CV_PI/180.)),2);
			double cg=pow(((i-400)-(ptg.x-400)*cos(angle*CV_PI/180.)-(400-ptg.y)*sin(angle*CV_PI/180.)),2);

			double bh=pow((sizeh.height*sin((angle-alphah)*CV_PI/180.)),2);
			double ah=pow((sizeh.width*cos((angle-alphah)*CV_PI/180.)),2);
			double ch=pow(((i-400)-(pth.x-400)*cos(angle*CV_PI/180.)-(400-pth.y)*sin(angle*CV_PI/180.)),2);

			double bi=pow((sizei.height*sin((angle-alphai)*CV_PI/180.)),2);
			double ai=pow((sizei.width*cos((angle-alphai)*CV_PI/180.)),2);
			double ci=pow(((i-400)-(pti.x-400)*cos(angle*CV_PI/180.)-(400-pti.y)*sin(angle*CV_PI/180.)),2);

			double bj=pow((sizej.height*sin((angle-alphaj)*CV_PI/180.)),2);
			double aj=pow((sizej.width*cos((angle-alphaj)*CV_PI/180.)),2);
			double cj=pow(((i-400)-(ptj.x-400)*cos(angle*CV_PI/180.)-(400-ptj.y)*sin(angle*CV_PI/180.)),2);

			//double t1=pa*2*sizea.height*sizea.width*sqrt(aa+ba-ca)/(aa+ba);
			double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
			if (aa+ba-ca<0) t1=0.;
			else t1=pa*2*sizea.height*sizea.width*sqrt(aa+ba-ca)/(aa+ba);

			if (ab+bb-cb<0) t2=0.;
			else t2=pb*2*sizeb.height*sizeb.width*sqrt(ab+bb-cb)/(ab+bb);
			
			if (ac+bc-cc<0) t3=0.;
			else t3=pc*2*sizec.height*sizec.width*sqrt(ac+bc-cc)/(ac+bc);
			
			if (ad+bd-cd<0) t4=0.;
			else t4=pd*2*sized.height*sized.width*sqrt(ad+bd-cd)/(ad+bd);

			if (ae+be-ce<0) t5=0.;
			else t5=pe*2*sizee.height*sizee.width*sqrt(ae+be-ce)/(ae+be);

			if (af+bf-cf<0) t6=0.;
			else t6=pf*2*sizef.height*sizef.width*sqrt(af+bf-cf)/(af+bf);
			
			if (ag+bg-cg<0) t7=0.;
			else t7=pg*2*sizeg.height*sizeg.width*sqrt(ag+bg-cg)/(ag+bg);
			
			if (ah+bh-ch<0) t8=0.;
			else t8=ph*2*sizeh.height*sizeh.width*sqrt(ah+bh-ch)/(ah+bh);

			if (ai+bi-ci<0) t9=0.;
			else t9=pi*2*sizei.height*sizei.width*sqrt(ai+bi-ci)/(ai+bi);

			if (aj+bj-cj<0) t10=0.;
			else t10=pj*2*sizej.height*sizej.width*sqrt(aj+bj-cj)/(aj+bj);
			 
			projfunc[i]=t1+t2+t3+t4+t5+t6+t7+t8+t9+t10;		//Now we got the back projection function!!!!!!!great!!!!

				/*+pb*2*sizeb.height*sizeb.width*sqrt(ab+bb-cb)/(ab+bb);
				/*+pc*2*sizec.height*sizec.width*sqrt(ac+bc-cc)/(ac+bc)
				+pd*2*sized.height*sized.width*sqrt(ad+bd-cd)/(ad+bd)
				+pe*2*sizee.height*sizee.width*sqrt(ae+be-ce)/(ae+be)
				+pf*2*sizef.height*sizef.width*sqrt(af+bf-cf)/(af+bf)
				+pg*2*sizeg.height*sizeg.width*sqrt(ag+bg-cg)/(ag+bg)
				+ph*2*sizeh.height*sizeh.width*sqrt(ah+bh-ch)/(ah+bh)
				+pi*2*sizei.height*sizei.width*sqrt(ai+bi-ci)/(ai+bi)
				+pj*2*sizej.height*sizej.width*sqrt(aj+bj-cj)/(aj+bj);*/
		}

		if (ctrl=='1')
		{
			if (ctrln=='y')	//there is noises