ctconv_reconstuct.cpp

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

// Plotting.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <cv.h>
#include <highgui.h>
#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <ctime>

using namespace std;

//the grey value 0~255 change to attenuation -4~2;
double grey2attenu(CvScalar newcolor)
{
	return (newcolor.val[0]/255*6-4);
}

//the attenuation -4~2 change to grey value 0~255;
CvScalar attenu2grey(double attenu)
{
	double a;
	a=(attenu+4)/6*255;
	return cvScalar(a);
}


int main()
{
	const char* saveimg1="S-L Head Model.bmp";
	const char* saveimg1rotate1="S-L Head Model Rotate.bmp";
	const char* saveimg1proj1="S-L Head Model Back Projection Reconstruction.bmp";

	const int imagedepth=IPL_DEPTH_8U;

	int SNR;		//the signal noise rate

	srand((unsigned)time( NULL ));		//set the random number seed system time
	
	int t,i,j,k;//loop indexes

	//the R-L convolution sequence (sharp frequence domain) 
	const unsigned int RLlen=800; 
	double rlconv[2*RLlen+1];
	for (i=0;i<2*RLlen+1;i++)
	{
		if (i==RLlen) rlconv[i]=0.25;
		if (((i%2)==0)&&(i!=RLlen)) rlconv[i]=0.;
		if (((i%2)!=0)&&(i!=RLlen)) rlconv[i]=-1.0/CV_PI/CV_PI/((i-RLlen)*(i-RLlen));
	}

	//the S-L convolution sequence (smooth frequence domain)
	const unsigned int SLlen=800;
	double slconv[2*SLlen+1];
	for (i=0;i<2*SLlen+1;i++) slconv[i]=-2./CV_PI/CV_PI/(4*(i-SLlen)*(i-SLlen)-1.);

	//my convolutions with different frequency respond
	const unsigned int MYCONVLEN=60;
	const double convsharphighpass[2*MYCONVLEN+1]={-0.000240858140043,-0.000133223717409,-0.000005154170445,0.000136939293637,0.000284263415617
,0.000424927462408,0.000543772703151,0.000623033280672,0.000643978280244,0.000589494970690
,0.000447359859918,0.000213744039480,-0.000103649544388,-0.000483483296816,-0.000890339445913
,-0.001276395732730,-0.001585422616489,-0.001758962282615,-0.001744177246757,-0.001502505408111
,-0.001017986023455,-0.000303978951053,0.000592974142356,0.001593185324658,0.002588399771881
,0.003451028970945,0.004047522961155,0.004254730417240,0.003977532969001,0.003165636034711
,0.001827224047075,0.000037291132652,-0.002061145017234,-0.004264068980025,-0.006321971804775
,-0.007963063485987,-0.008922404125941,-0.008974258625376,-0.007964308634454,-0.005838020030214
,-0.002661532295388,0.001368067994971,0.005925329190586,0.010572886030539,0.014791359623433
,0.018020790955999,0.019710305412978,0.019371587456048,0.016631035912732,0.011275272117329
,0.003285033288548,-0.007146610768674,-0.019614402582766,-0.033520369782713,-0.048111703314857
,-0.062534256227342,-0.075897240654455,-0.087343467603161,-0.096118785174237,-0.101634309392224
,0.896484375000000,-0.101634309392238,-0.096118785174230,-0.087343467603166,-0.075897240654451
,-0.062534256227345,-0.048111703314854,-0.033520369782715,-0.019614402582764,-0.007146610768675
,0.003285033288550,0.011275272117327,0.016631035912733,0.019371587456048,0.019710305412979
,0.018020790955998,0.014791359623433,0.010572886030538,0.005925329190586,0.001368067994970
,-0.002661532295387,-0.005838020030214,-0.007964308634453,-0.008974258625377,-0.008922404125941
,-0.007963063485987,-0.006321971804775,-0.004264068980025,-0.002061145017234,0.000037291132652
,0.001827224047075,0.003165636034711,0.003977532969001,0.004254730417240,0.004047522961155
,0.003451028970945,0.002588399771881,0.001593185324658,0.000592974142356,-0.000303978951053
,-0.001017986023455,-0.001502505408111,-0.001744177246757,-0.001758962282615,-0.001585422616488
,-0.001276395732730,-0.000890339445913,-0.000483483296816,-0.000103649544388,0.000213744039480
,0.000447359859918,0.000589494970690,0.000643978280244,0.000623033280672,0.000543772703151
,0.000424927462408,0.000284263415617,0.000136939293637,-0.000005154170445,-0.000133223717409
,-0.000240858140043};

	const double convsmoothhighpass[2*MYCONVLEN+1]={-0.000002427537616,-0.000012715444644,-0.000019791050885,-0.000011820106991,0.000011547521970
,0.000035523379635,0.000040431663449,0.000019333465792,-0.000012824054095,-0.000030254685927
,-0.000021089278352,-0.000002283973376,-0.000005132948939,-0.000039619622270,-0.000073044237627
,-0.000052991765392,0.000037428432917,0.000145402673537,0.000180863925411,0.000098914398005
,-0.000045175199632,-0.000135268947737,-0.000106178639084,-0.000019850593455,-0.000012239662765
,-0.000141269406895,-0.000286808994387,-0.000233350508991,0.000107290068507,0.000545731825186
,0.000725137770416,0.000440211849515,-0.000127537871162,-0.000528610804361,-0.000465660928969
,-0.000122471730292,-0.000021075810004,-0.000459221876739,-0.001060321181827,-0.000980908486490
,0.000263438558135,0.002088577520400,0.003065308588043,0.002112655141408,-0.000333635051571
,-0.002421651134172,-0.002490909020635,-0.000881487836014,-0.000028254265008,-0.002150455179072
,-0.006191072210049,-0.007079981439436,0.000802101130858,0.017521089286055,0.033508332582753
,0.032000992750966,-0.001608723022364,-0.068619506600821,-0.151808166367031,-0.221087162237070
,0.751953125000000,-0.221087162237083,-0.151808166367023,-0.068619506600825,-0.001608723022359
,0.032000992750964,0.033508332582755,0.017521089286053,0.000802101130860,-0.007079981439438
,-0.006191072210047,-0.002150455179073,-0.000028254265007,-0.000881487836015,-0.002490909020634
,-0.002421651134173,-0.000333635051570,0.002112655141407,0.003065308588044,0.002088577520400
,0.000263438558136,-0.000980908486490,-0.001060321181826,-0.000459221876739,-0.000021075810004
,-0.000122471730292,-0.000465660928969,-0.000528610804361,-0.000127537871162,0.000440211849515
,0.000725137770417,0.000545731825185,0.000107290068507,-0.000233350508991,-0.000286808994386
,-0.000141269406895,-0.000012239662765,-0.000019850593455,-0.000106178639084,-0.000135268947737
,-0.000045175199632,0.000098914398005,0.000180863925412,0.000145402673537,0.000037428432917
,-0.000052991765392,-0.000073044237627,-0.000039619622270,-0.000005132948939,-0.000002283973376
,-0.000021089278352,-0.000030254685927,-0.000012824054095,0.000019333465792,0.000040431663449
,0.000035523379635,0.000011547521970,-0.000011820106991,-0.000019791050885,-0.000012715444644
,-0.000002427537616};

	const double convalteredhighpass[2*MYCONVLEN+1]={0.000000181027643,-0.000003942403032,0.000006876573808,-0.000001311943457,0.000001240927289
,-0.000015630412142,-0.000008084247087,-0.000013639044966,0.000009579609037,0.000001006716889
,0.000012792667047,-0.000016120244012,0.000000382204222,-0.000014740213506,0.000025976037004
,-0.000003936750831,0.000007514361164,-0.000065227661497,-0.000036171867164,-0.000062972386613
,0.000038420754771,0.000003832841535,0.000057331448321,-0.000066258957058,0.000000910336743
,-0.000064239761697,0.000105166307622,-0.000010997232752,0.000039539646523,-0.000251018348747
,-0.000145032599641,-0.000257828991981,0.000134840622810,0.000011506959192,0.000228306846623
,-0.000245194047669,0.000001566274690,-0.000260837700374,0.000404046795343,-0.000025825522444
,0.000196471604766,-0.000991246483227,-0.000613076193390,-0.001155444673814,0.000532397013646
,0.000032493838918,0.001126601521517,-0.001185479786485,0.000002098748393,-0.001566879254362
,0.002472823042720,-0.000070895528478,0.001827449143686,-0.008641767562886,-0.006701688147097
,-0.016539142045600,0.008748941427845,0.000204846139101,0.064168462375733,-0.213831309596134
,0.350585937500000,-0.213831309596139,0.064168462375740,0.000204846139097,0.008748941427849
,-0.016539142045603,-0.006701688147095,-0.008641767562888,0.001827449143688,-0.000070895528480
,0.002472823042721,-0.001566879254363,0.000002098748394,-0.001185479786486,0.001126601521518
,0.000032493838917,0.000532397013647,-0.001155444673815,-0.000613076193390,-0.000991246483228
,0.000196471604766,-0.000025825522444,0.000404046795344,-0.000260837700374,0.000001566274690
,-0.000245194047669,0.000228306846624,0.000011506959192,0.000134840622810,-0.000257828991981
,-0.000145032599641,-0.000251018348747,0.000039539646523,-0.000010997232752,0.000105166307622
,-0.000064239761697,0.000000910336743,-0.000066258957058,0.000057331448322,0.000003832841535
,0.000038420754771,-0.000062972386613,-0.000036171867164,-0.000065227661497,0.000007514361164
,-0.000003936750831,0.000025976037004,-0.000014740213506,0.000000382204222,-0.000016120244012
,0.000012792667047,0.000001006716889,0.000009579609037,-0.000013639044966,-0.000008084247087
,-0.000015630412142,0.000001240927289,-0.000001311943457,0.000006876573808,-0.000003942403032
,0.000000181027643};
	
	/*{-0.000029763026467,0.000000192015897,-0.000034952971794,0.000062406360220,-0.000004721311561
,0.000042828815164,-0.000256640222410,-0.000186699770128,-0.000408861849578,0.000216112176921
,0.000014940721268,0.000579757223187,-0.000675127348337,0.000001308936557,-0.001059823739309
,0.001797614811989,-0.000054922736591,0.001496778308186,-0.007427379416449,-0.006001271632594
-0.015326369674793,0.008334919521694,0.000199358069036,0.063400881729293,-0.213190090694259
,0.350585937500000,-0.213190090694259,0.063400881729293,0.000199358069036,0.008334919521694
,-0.015326369674793,-0.006001271632594,-0.007427379416449,0.001496778308186,-0.000054922736591
,0.001797614811989,-0.001059823739309,0.000001308936557,-0.000675127348337,0.000579757223187
,0.000014940721268,0.000216112176921,-0.000408861849578,-0.000186699770128,-0.000256640222410
,0.000042828815164,-0.000004721311561,0.000062406360220,-0.000034952971794,0.000000192015897
,-0.000029763026467};*/ 

	/*{-0.000027589289829,0.000000363573931,-0.000028008834174,0.000076127195701,-0.000013816194114,
-0.000053474877605,-0.000500308294097,-0.000560726047498,-0.000782923521239,0.000013430148016,
0.000045832378076,0.000720432623783,-0.000606790719984,0.000002482521439,-0.000885560227526,
0.002215263889149,-0.000162900209098,-0.000996210633479,-0.014416121490614,-0.018006036765922,
-0.029571872507624,-0.002071553410706,0.000600530429604,0.078612035788222,-0.185093350735735,
0.383789062500000,-0.185093350735735,0.078612035788222,0.000600530429604,-0.002071553410706,
-0.029571872507624,-0.018006036765922,-0.014416121490614,-0.000996210633479,-0.000162900209098,
0.002215263889149,-0.000885560227526,0.000002482521439,-0.000606790719984,0.000720432623783,
0.000045832378076,0.000013430148016,-0.000782923521239,-0.000560726047498,-0.000500308294097,
-0.000053474877605,-0.000013816194114,0.000076127195701,-0.000028008834174,0.000000363573931,
-0.000027589289829};*/ 


	char ctrl;	//can give 4 method: 
				//ctrl=0-->direct back projection 
				//ctrl=1-->R-L convolution 
				//ctrl=2-->S-L convolution 
				//ctrl=3-->my convolution 
	char ctrln;	//1-->with noises
					//0-->no noises

	char ctrlmy;	//ctrl=0-->sharp high pass convolution
				//ctrl=1-->smooth high pass convolution 
				//ctrl=2-->altered high pass convolution 

	char ctrlnew;	//this controls whether to enter the parameters of new objects

	char ctrlquit;	//this controls whether to quit the programme

	double totaltime;	//calculate the reconstructing time


	do{
	//------------------------------------------//
	cout<<"---Welcome using this simple CT reconstruction simulation software---"<<endl;
	cout<<"-----------------------Version 1.00 By Simon-------------------------"<<endl<<endl;
	cout<<"Show the original image"<<endl;
	
	//define the size of the image
	const int width=800;
	const int height=800;
	const int N=height*width;
	
	//define the background grey value 
	const int Background=0; //0~255
	
	//define some center points to plot the model
	//------------------------------------------//
	CvPoint pta;
	pta.x=400;
	pta.y=400;

	CvPoint ptb;
	ptb.x=400;
	ptb.y=407.36;

	CvPoint ptc;
	ptc.x=488;
	ptc.y=400;

	CvPoint ptd;
	ptd.x=312;
	ptd.y=400;

	CvPoint pte;
	pte.x=400;
	pte.y=260;

	CvPoint ptf;
	ptf.x=400;
	ptf.y=360;

	CvPoint ptg;
	ptg.x=400;
	ptg.y=440;

	CvPoint pth;
	pth.x=368;
	pth.y=642;

	CvPoint pti;
	pti.x=400;
	pti.y=642;

	CvPoint ptj;
	ptj.x=424;
	ptj.y=642;
	//------------------------------------------//

	//define some grey value to describe the attenuation of the tissue
	//------------------------------------------//
	CvScalar colora=attenu2grey(2.0);	
	CvScalar colorb=attenu2grey(-0.98);
	CvScalar colorc=attenu2grey(-0.02);
	CvScalar colord=attenu2grey(-0.02);
	CvScalar colore=attenu2grey(0.01);
	CvScalar colorf=attenu2grey(0.01);
	CvScalar colorg=attenu2grey(0.02);
	CvScalar colorh=attenu2grey(0.01);
	CvScalar colori=attenu2grey(0.01);
	CvScalar colorj=attenu2grey(0.01);
	//------------------------------------------//

	//define the size of the ellipse objects
	//------------------------------------------//
	CvSize sizea=cvSize(368,276);
	CvSize sizeb=cvSize(349.6,264.96);
	CvSize sizec=cvSize(124,44);
	CvSize sized=cvSize(164,64);
	CvSize sizee=cvSize(100,84);
	CvSize sizef=cvSize(18.4,18.4);
	CvSize sizeg=cvSize(18.4,18.4);
	CvSize sizeh=cvSize(18.4,9.2);
	CvSize sizei=cvSize(9.2,9.2);
	CvSize sizej=cvSize(18.4,9.2);

	//define the attenuation of the ellipse objects	
	//unmodified model
	//double pa=2.,pb=-0.98,pc=-0.01,pd=-0.02,pe=0.01,pf=0.01,pg=0.02,ph=0.01,pi=0.01,pj=0.01;

	//modified model
	double pa=1.,pb=-0.98,pc=-0.1,pd=-0.2,pe=0.1,pf=0.1,pg=0.2,ph=0.1,pi=0.1,pj=0.1;