Rev 1 | Blame | Compare with Previous | Last modification | View Log | RSS feed
//// ../bin/FBP2D FBP2D.par//// ********************************************************************// * License and Disclaimer *// * *// * The GAMOS software is copyright of the Copyright Holders of *// * the GAMOS Collaboration. It is provided under the terms and *// * conditions of the GAMOS Software License, included in the file *// * LICENSE and available at http://fismed.ciemat.es/GAMOS/license .*// * These include a list of copyright holders. *// * *// * Neither the authors of this software system, nor their employing *// * institutes,nor the agencies providing financial support for this *// * work make any representation or warranty, express or implied, *// * regarding this software system or assume any liability for its *// * use. Please see the license in the file LICENSE and URL above *// * for the full disclaimer and the limitation of liability. *// * *// * This code implementation is the result of the scientific and *// * technical work of the GAMOS collaboration. *// * By using, copying, modifying or distributing the software (or *// * any work based on the software) you agree to acknowledge its *// * use in resulting scientific publications, and indicate your *// * acceptance of all terms of the GAMOS Software license. *// ********************************************************************//#include "PETProjDataMgr.h"#include "TH2F.h"#include "TH3F.h"#include "TMath.h"#include "TRandom3.h"#include <iomanip>#include <iostream>#include <cstdlib>using namespace std;//----------------------------------------------------------------------PETProjDataMgr* PETProjDataMgr::m_Instance = 0;//----------------------------------------------------------------------PETProjDataMgr* PETProjDataMgr::GetInstance(){if( !m_Instance ) {m_Instance = new PETProjDataMgr;}return m_Instance;}//-----------------------------------------------------------------------PETProjDataMgr::PETProjDataMgr(){/*// ARGUMENTS:" cout << " -------------------------- \n"" Arguments convention: \n"" -a Axial FOV (mm), <theDist_axial=100.0> \n"" -d Diameter Transaxial FOV (mm), <m_RingDiameter=300.0> \n"" -i Type of the input file (by default: 0 = Arce_binary), <typeINfile=0> \n"" \n"" -m Maximum ring difference (by default: -1 = m_NOfPlanes), <m_MaxRingDifferenceiff> \n"" -n Name of output file, <m_Filename> \n"" -p Axial number of planes, <m_NOfPlanes> \n"" -r Number of bins, \"distancias\", <m_NOfBins> \n"// -s Span TO DO:span !!!!!" -t Number of angular views, \"direcciones\", <m_NOfAngles> \n"" -v Verbosity (by default: 0=silent, 3=debug), <verbos> \n"" -x Maximum number of coincidences to be stored (by default: -1 = no limit), <Max_Coinci> \n"" -o Output type (by default: 0 = mcc Interfile, 1 = STIR Interfile), <OutType> \n"" \n"" PET Reconstruction. CIEMAT 2009-11 \n"" mario.canadas@ciemat.es \n"" -------------------------- \n";*/m_rnd= new TRandom3();m_AxialDistance = (9-1)*2.25; // Axial pixel dimension*NOfPlanesm_RingDiameter = 120.0; // notranji premer petam_NOfPlanes = 9; // stevilo ravninm_NOfBins = 128; // stevilo binov v razdaljim_nch = 128; // stevilo padov okoli in okolim_NOfAngles = TMath::Abs(m_nch)/2; // stevilo kotov = stevilo padov okoli in okoli /2m_MaxRingDifference = -1;//m_MaxRingDifference = 3; // najvecja razdalja med padi// toDo: theSpan = int(GmParameterMgr::GetInstance()->GetNumericValue("PET:ProjData:Span",1));m_OutFormat = 1; // 1.. projections 0 .. imagem_Filename = "sino3D";m_Debug=1;if (m_MaxRingDifference==-1) m_MaxRingDifference=m_NOfPlanes-1;m_TotalAxialPlanes=m_NOfPlanes*m_NOfPlanes;if (m_OutFormat==1) m_TotalAxialPlanes= (2*m_NOfPlanes-1 - m_MaxRingDifference)*m_MaxRingDifference + m_NOfPlanes; // total number of Axial planes (segments*planes) in STIR format/*--- Initialize sino3D ---*/m_projections = new SINO_TYPE**[m_NOfBins];for(int i=0;i<m_NOfBins;i++){m_projections[i] = new SINO_TYPE*[m_NOfAngles];for(int j=0;j<m_NOfAngles;j++){m_projections[i][j] = new SINO_TYPE[m_TotalAxialPlanes]; /// ! If m_OutFormat==1 (STIR output):Matrix size depends on the MAX_Ring_Differencefor(int k=0;k<m_TotalAxialPlanes;k++){m_projections[i][j][k]=0;}}}m_TotalProjectionCoincidences=0;m_TotalCoincidences=0;//OutputType = "pet";}//-----------------------------------------------------------------------void PETProjDataMgr::SetProjection( int axialplane, TH2F * h){for(int i=0;i<m_NOfBins;i++){for(int j=0;j<m_NOfAngles;j++){m_projections[i][j][axialplane]=h->GetBinContent(i+1,j+1);}}}//-----------------------------------------------------------------------//-----------------------------------------------------------// from Gate//------------------------------------------------------------double ComputeSinogramS(double X1, double Y1, double X2, double Y2){double s;double denom = (Y1-Y2) * (Y1-Y2) + (X2-X1) * (X2-X1);if (denom!=0.) {denom = sqrt(denom);s = ( X1 * (Y2-Y1) + Y1 * (X1-X2) ) / denom;} else {s = 0.;}double theta;if ((X1-X2)!=0.) {theta=atan((X1-X2) /(Y1-Y2));} else {theta=3.1416/2.;}if ((theta > 0.) && ((X1-X2) > 0.)) s = -s;if ((theta < 0.) && ((X1-X2) < 0.)) s = -s;if ( theta < 0.) {theta = theta+3.1416;s = -s;}return s;}void PETProjDataMgr::AddEvent( const TVector3& pos1 , const TVector3& pos2){int z1_i, z2_i;//for discretization on the crystal: int x1_i, x2_i, y1_i, y2_i;double z1_abs=pos1.z()+m_AxialDistance/2;double z2_abs=pos2.z()+m_AxialDistance/2;double a, b, phi, dis;int phi_i, dis_i;int ring_diff;double _PI=2*asin(1);m_TotalCoincidences++;z1_i=(int)(m_NOfPlanes* z1_abs/m_AxialDistance); //round --> m_NOfPlanes+1 ...z2_i=(int)(m_NOfPlanes* z2_abs/m_AxialDistance);// control; if z_i out of range: returnif ( (pos1.x()==pos2.x()) && (pos1.y()==pos2.y()) ) {#ifndef GAMOS_NO_VERBOSEif( m_Debug ) {cout << "PETProjDataMgr::AddEvent:WARNING! Event_1 == Event_2 ; x= " << pos2.x() << ", y= " << pos2.y() << ", z= " << pos2.z() << endl;}#endifreturn;}if ( (z1_i<0) || (z2_i<0) || (z1_i>= m_NOfPlanes) || (z2_i>= m_NOfPlanes) ) {#ifndef GAMOS_NO_VERBOSEif( m_Debug ) {cout << "PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Axial): x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl;}#endifreturn;}ring_diff = (int)fabs(z1_i-z2_i);// max ring difference; control:if (ring_diff > m_MaxRingDifference) {#ifndef GAMOS_NO_VERBOSEif( m_Debug ) {cout <<"PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Max. Ring Diff.): " << ring_diff << ">" << m_MaxRingDifference << " x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl;}#endifreturn;}a=(double)(pos2.y()- pos1.y());b=(double)(pos2.x()- pos1.x());if (a==0.0){phi=_PI*0.5;}else{phi=atan(b/a);}if (phi<0) phi = phi +_PI;dis=pos1.x()*cos(phi) - pos1.y()*sin(phi);//dis=ComputeSinogramS(pos1.x(), pos1.y(), pos2.x(), pos2.x());// control; transaxial FOVif ( fabs(dis) > m_RingDiameter*0.5 ) {#ifndef GAMOS_NO_VERBOSEif( m_Debug ) {cout << "PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Transaxial): x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl;}#endifreturn;}dis = dis + m_RingDiameter*0.5;// discret values:phi_i=(int)round( (double)(m_NOfAngles-1)*phi/_PI );dis_i=(int)round( (double)(m_NOfBins-1)*dis/(double)m_RingDiameter );if ((phi_i>=m_NOfAngles) || (dis_i>=m_NOfBins)) return; // only possible "=" because 'round' check it..// OLD: (SRRB included) sino3D[dis_i][phi_i][ (z1_i+z2_i)+ring_diff*(m_NOfPlanes-1) ]++;int Zpos;if (m_OutFormat==0) {Zpos = (z1_i*m_NOfPlanes + z2_i);}else{if (z1_i>=z2_i) { // SIN Max Ring_Diff: Zpos= ( ((m_NOfPlanes-ring_diff)*(m_NOfPlanes-1-ring_diff))/2 + z2_i );Zpos= ( ((2*m_NOfPlanes-1 - m_MaxRingDifference - ring_diff)*(m_MaxRingDifference - ring_diff))/2 + z2_i);}else{Zpos= ( (m_TotalAxialPlanes) - ((2*m_NOfPlanes-1 - m_MaxRingDifference - ring_diff +1)*(m_MaxRingDifference - ring_diff +1))/2 + z1_i );}}m_projections[dis_i][phi_i][ Zpos ]++;m_TotalProjectionCoincidences++;#ifndef GAMOS_NO_VERBOSEif( m_Debug >1) {cout << "PETProjDataMgr::AddEvent: x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl;cout << "PETProjDataMgr::AddEvent: Sinogram pos.: distance(s)= " << dis_i << ", angular view(phi)= " << phi_i << " ; Zpos=" << Zpos <<"; Segment (Ring diff.) = " << ring_diff << endl;}#endif}//-----------------------------------------------------------------------PETProjDataMgr::~PETProjDataMgr(){int i,j;for(i=0;i<m_NOfBins;i++){for(j=0;j<m_NOfAngles;j++){free(m_projections[i][j]);}free(m_projections[i]);}free(m_projections);}/* TO DO: call lm_to_sino3D program//-----------------------------------------------------------------------void PETIOMgr::ReadFile(){if( !theFileIn ) OpenFileIn();PETOutput po;G4bool bEof;for(;;) {po = ReadEvent( bEof );// theFileIn->read(reinterpret_cast<char *>(&po),sizeof(PetOutput));if(bDumpCout) PrintEvent(" PETOutput: ", po, bCartesian);if( bEof ) break;}}//-----------------------------------------------------------------------PETOutput PETProjDataMgr::ReadEvent( G4bool& bEof ){if( theFileIn == 0 ){G4Exception("PETIOMgr::ReadEvent, file not opened, call OpenFileIn() first ");}PETOutput po;fread (&po, sizeof(struct PETOutput),1,theFileIn);if ( feof( theFileIn ) ) {bEof = TRUE;} else {bEof = FALSE;}return po;}*///-----------------------------------------------------------------------void PETProjDataMgr::WriteInterfile(){char name_hv[512];char name_v[512];if (m_OutFormat==0){strcpy(name_hv, m_Filename);strcpy(name_v,m_Filename);strcat(name_hv, ".hv");strcat(name_v, ".v");fp=fopen(name_hv, "w");fprintf (fp, "!INTERFILE := \n");fprintf (fp, "name of data file := %s\n", name_v);fprintf (fp, "!GENERAL DATA := \n");fprintf (fp, "!GENERAL IMAGE DATA :=\n");fprintf (fp, "!type of data := tomographic\n");fprintf (fp, "!version of keys := 3.3\n");fprintf (fp, "!data offset in bytes := 0\n");fprintf (fp, "imagedata byte order := littleendian\n");fprintf (fp, "!PET STUDY (General) :=\n");fprintf (fp, "!PET data type := 3D-Sinogram\n");fprintf (fp, "process status := Reconstructed\n");fprintf (fp, "!number format := unsigned short\n");fprintf (fp, "!number of bytes per pixel := 2\n");fprintf (fp, "number of dimensions := 3\n");fprintf (fp, "matrix axis label [1] := x\n");fprintf (fp, "!matrix size [1] := %i\n",m_NOfBins);fprintf (fp, "scaling factor (mm/pixel) [1] := %f\n",(float)(m_RingDiameter/(m_NOfBins-1.0)));fprintf (fp, "matrix axis label [2] := y\n");fprintf (fp, "!matrix size [2] := %i\n",m_NOfAngles);fprintf (fp, "scaling factor (degree/pixel) [2] := %f\n",(float)(360./m_NOfAngles));fprintf (fp, "matrix axis label [3] := z\n");fprintf (fp, "!matrix size [3] := %i\n",m_NOfPlanes*m_NOfPlanes);fprintf (fp, "scaling factor (mm/pixel) [3] := %f\n",(float)(m_AxialDistance/(m_NOfPlanes-1.0)));fprintf (fp, "number of slices := %i\n",m_NOfPlanes*m_NOfPlanes);fprintf (fp, "number of time frames := 1\n");fprintf (fp, "image scaling factor[1] := 1\n");fprintf (fp, "data offset in bytes[1] := 0\n");fprintf (fp, "quantification units := 1\n");fprintf (fp, "!END OF INTERFILE := \n");fclose(fp);//(size_t)(m_NOfBins*m_NOfAngles*m_NOfPlanes*m_NOfPlanes);}else{strcpy(name_hv, m_Filename);strcpy(name_v,m_Filename);strcat(name_hv, ".hs"); // STIR extension: .hs .sstrcat(name_v, ".s");fp=fopen(name_hv, "w");fprintf (fp, "!INTERFILE := \n");fprintf (fp, "name of data file := %s\n",name_v);fprintf (fp, "!GENERAL DATA := \n");fprintf (fp, "!GENERAL IMAGE DATA :=\n");fprintf (fp, "!type of data := PET\n");// fprintf (fp, "!version of keys := 3.3\n"); STIR format is not 3.3 (almost but not completely), ERROR in STIR if it is not removed// fprintf (fp, "!data offset in bytes := 0\n");fprintf (fp, "imagedata byte order := littleendian\n");fprintf (fp, "!PET STUDY (General) :=\n");fprintf (fp, "!PET data type := Emission\n");fprintf (fp, "applied corrections := {arc correction}\n"); // {none}\n");// fprintf (fp, "process status := Reconstructed\n");fprintf (fp, "!number format := unsigned integer\n");fprintf (fp, "!number of bytes per pixel := 2\n");fprintf (fp, "number of dimensions := 4\n");fprintf (fp, "matrix axis label [4] := segment\n");fprintf (fp, "!matrix size [4] := %i\n",m_MaxRingDifference*2 + 1);// fprintf (fp, "scaling factor (mm/pixel) [1] := %f\n",(float)(d_FOV/(m_NOfBins-1)));fprintf (fp, "matrix axis label [3] := axial coordinate\n");fprintf (fp, "!matrix size [3] := { ");if (m_MaxRingDifference==0){fprintf (fp, "%i}\n", m_NOfPlanes);}else{for(int m=m_NOfPlanes-m_MaxRingDifference;m<=m_NOfPlanes;m++) fprintf (fp, "%i,", m);for(int m=m_NOfPlanes-1;m>m_NOfPlanes-m_MaxRingDifference;m--) fprintf (fp, "%i,", m);fprintf (fp, "%i}\n", m_NOfPlanes-m_MaxRingDifference);}fprintf (fp, "matrix axis label [2] := view\n");fprintf (fp, "!matrix size [2] := %i\n",m_NOfAngles);fprintf (fp, "matrix axis label [1] := tangential coordinate\n");fprintf (fp, "!matrix size [1] := %i\n",m_NOfBins);fprintf (fp, "minimum ring difference per segment := {"); // TO DO : add SPAN (m_MaxRingDifferenceiff per seg. variable)fprintf (fp, "%i", -m_MaxRingDifference);for(int m=-m_MaxRingDifference+1;m<=m_MaxRingDifference;m++) fprintf (fp, ",%i", m);fprintf (fp, "}\n");fprintf (fp, "maximum ring difference per segment := {"); // TO DO : add SPAN (m_MaxRingDifferenceiff per seg. variable)fprintf (fp, "%i", -m_MaxRingDifference);for(int m=-m_MaxRingDifference+1;m<=m_MaxRingDifference;m++) fprintf (fp, ",%i", m);fprintf (fp, "}\n");fprintf (fp, "inner ring diameter (cm) := %f\n", m_RingDiameter/10); // STIR Required parameter, now assigned to FOV (not detectors)fprintf (fp, "average depth of interaction (cm) := 0.0001\n");fprintf (fp, "default bin size (cm) := %f\n",0.1*((float)m_RingDiameter/((float)m_NOfBins-1.0)) );fprintf (fp, "number of rings := %i\n",m_NOfPlanes );fprintf (fp, "distance between rings (cm) := %f\n", 0.1*((float)m_AxialDistance/(float)(m_NOfPlanes-1)) ); // Axial pixel dimensionfprintf (fp, "number of detectors per ring := %i\n",m_NOfAngles*2 );// fprintf (fp, "number of slices := %i\n",m_NOfPlanes*m_NOfPlanes);fprintf (fp, "number of time frames := 1\n");fprintf (fp, "image scaling factor[1] := 1\n");fprintf (fp, "data offset in bytes[1] := 0\n");fprintf (fp, "quantification units := 1\n");fprintf (fp, "!END OF INTERFILE := \n");fclose(fp);}m_Buffer = (SINO_TYPE*) malloc( m_NOfBins*m_NOfAngles*m_TotalAxialPlanes*sizeof(SINO_TYPE));long unsigned int cont=0;int i,j,k;for(k=0;k<m_TotalAxialPlanes;k++){for(j=0;j<m_NOfAngles;j++){for(i=0;i<m_NOfBins;i++){m_Buffer[cont]=m_projections[i][j][k];cont++;}}}fp=fopen(name_v, "wb");//cout << 4096*sizeof(SINO_TYPE) << endl;int nb=fwrite(m_Buffer,1,m_NOfBins*m_NOfAngles*m_TotalAxialPlanes*sizeof(SINO_TYPE), fp);fclose(fp);#ifndef GAMOS_NO_VERBOSEcout << "PETProjDataMgr::WriteInterfile: File name: " << m_Filename << endl;cout << "PETProjDataMgr::WriteInterfile: Numer of bytes written: " << nb << endl;cout << "PETProjDataMgr::WriteInterfile: Planes = " << m_NOfPlanes << "; bins = " << m_NOfBins << "; ang_views = " << m_NOfAngles << endl;cout << "PETProjDataMgr::WriteInterfile: Dimensions (mm): Transaxial FOV = " << m_RingDiameter << "; Axial FOV = " << m_AxialDistance << " ; Transaxial_pix = " << m_RingDiameter/(m_NOfBins-1) <<"; Plane width = " << m_AxialDistance/(m_NOfPlanes-1) << endl; // Image Axial Pixel(ssrb) == 0.5*(Plane_Width);cout << "... " << endl;cout << "PETProjDataMgr::WriteInterfile: Total Coinci: " << m_TotalCoincidences << endl;cout << "PETProjDataMgr::WriteInterfile: Sino3D Coinci: " << m_TotalProjectionCoincidences << endl;#endif}TVector3 PETProjDataMgr::Hits2Digits(const TVector3 &r){if (!m_nch) return r;float smear=0.5;if (m_nch<0) smear=m_rnd->Rndm();double angle = TMath::ATan2(r.X(),r.Y()); // vrne kot med -pi in piif (angle<0) angle+=TMath::TwoPi();angle= (int(angle/TMath::TwoPi()*TMath::Abs(m_nch))+smear)*TMath::TwoPi()/TMath::Abs(m_nch);//(m_rnd->Rndm()-0.5)*m_AxialDistance;return TVector3(sin(angle), cos(angle),0); // z coordinata ni cisto v redu}int PETProjDataMgr::FwdProject(double x,double y, double z, int nmax, TH1*h){TVector3 r(x,y,z);int h2d=h->InheritsFrom("TH2F");double tfac=m_RingDiameter*m_RingDiameter/4-r.Mag2();double rfac= m_AxialDistance/m_RingDiameter;for (int i=0;i<nmax;i++){double phi= m_rnd->Rndm()*TMath::Pi();TVector3 s(1,0,0);s.SetPhi(phi);double sign = (m_rnd->Rndm()>0.5)? 1 : 0;double theta = TMath::ACos(m_rnd->Rndm()*rfac);theta+=sign*TMath::Pi();s.SetTheta(theta);double t=r*s;TVector3 rx=r-t*s;double d=TMath::Sqrt(t*t+tfac);TVector3 r1=rx+d*s;TVector3 r2=rx-d*s;//r1=Hits2Digits(r1);//r2=Hits2Digits(r2);AddEvent( r1 , r2);if (h!=NULL){TVector3 s1=r2-r1;double s1len= s1.Mag();int niter=int (100*s1len/m_RingDiameter);for (int j=0;j<niter;j++){r2=r1+m_rnd->Rndm()*s1;if (h2d) ((TH2F *) h)->Fill(r2.X(),r2.Y());else ((TH3F *) h)->Fill(r2.X(),r2.Y(),r2.Z());}}}return 0;}int PETProjDataMgr::FwdProject(TH2F *img, TH2F *h){for (int i=0;i<img->GetNbinsX();i++) {double x_=img->GetXaxis()->GetBinCenter( i+1 );for (int j=0;j<img->GetNbinsY();j++) {double y_=img->GetYaxis()->GetBinCenter( j+1 );double density= img->GetBinContent(i+1,j+1);if (density>0) FwdProject(x_,y_,m_AxialDistance*(m_rnd->Rndm()-0.5), density,h);}}return 0;}int PETProjDataMgr::FwdProject(TH3F *img, TH3F *h){for (int i=0;i<img->GetNbinsX();i++) {double x_=img->GetXaxis()->GetBinCenter( i+1 );for (int j=0;j<img->GetNbinsY();j++) {double y_=img->GetYaxis()->GetBinCenter( j+1 );for (int k=0;k<img->GetNbinsZ();k++) {double z_=img->GetZaxis()->GetBinCenter( k+1 );double density= img->GetBinContent(i+1,j+1,k+1);if (density>0) FwdProject(x_,y_,z_, density,h);}}}return 0;}TH2F *PETProjDataMgr::Phantom(int kaj){TH2F *img= new TH2F("img","Original Image",100,-50,50,100,-50,50);// izberi sliko 0: kroglice, 1: point source 2: central ballswitch (kaj){case 0:for (int i=0;i<img->GetNbinsX();i++) {for (int j=0;j<img->GetNbinsY();j++) {double x_=img->GetXaxis()->GetBinCenter( i+1 );double y_=img->GetYaxis()->GetBinCenter( j+1 );double density=1000;if ((x_*x_+y_*y_)<6) img->SetBinContent(i+1,j+1,density);density=500; if ((x_-25)*(x_-25)+y_*y_<12) img->SetBinContent(i+1,j+1,density);density=2000; if ((y_-25)*(y_-25)+x_*x_<2) img->SetBinContent(i+1,j+1,density);}}break;case 2:for (int i=0;i<img->GetNbinsX();i++) {for (int j=0;j<img->GetNbinsY();j++) {double x_=img->GetXaxis()->GetBinCenter( i+1 );double y_=img->GetYaxis()->GetBinCenter( j+1 );double density=1000;if ((x_*x_+y_*y_)<12.5) img->SetBinContent(i+1,j+1,density);}}break;case 1:img->Fill(25,25,10000);break;}return img;}