Subversion Repositories f9daq

Rev

Rev 1 | Blame | Compare with Previous | Last modification | View Log | RSS feed

  1. //
  2. // ../bin/FBP2D FBP2D.par
  3. //
  4. // ********************************************************************
  5. // * License and Disclaimer                                           *
  6. // *                                                                  *
  7. // * The  GAMOS software  is  copyright of the Copyright  Holders  of *
  8. // * the GAMOS Collaboration.  It is provided  under  the  terms  and *
  9. // * conditions of the GAMOS Software License,  included in the  file *
  10. // * LICENSE and available at  http://fismed.ciemat.es/GAMOS/license .*
  11. // * These include a list of copyright holders.                       *
  12. // *                                                                  *
  13. // * Neither the authors of this software system, nor their employing *
  14. // * institutes,nor the agencies providing financial support for this *
  15. // * work  make  any representation or  warranty, express or implied, *
  16. // * regarding  this  software system or assume any liability for its *
  17. // * use.  Please see the license in the file  LICENSE  and URL above *
  18. // * for the full disclaimer and the limitation of liability.         *
  19. // *                                                                  *
  20. // * This  code  implementation is the result of  the  scientific and *
  21. // * technical work of the GAMOS collaboration.                       *
  22. // * By using,  copying,  modifying or  distributing the software (or *
  23. // * any work based  on the software)  you  agree  to acknowledge its *
  24. // * use  in  resulting  scientific  publications,  and indicate your *
  25. // * acceptance of all terms of the GAMOS Software license.           *
  26. // ********************************************************************
  27. //
  28. #include "PETProjDataMgr.h"
  29. #include "TH2F.h"
  30. #include "TH3F.h"
  31. #include "TMath.h"
  32. #include "TRandom3.h"
  33.  
  34. #include <iomanip>
  35. #include <iostream>
  36. #include <cstdlib>
  37.  
  38. using namespace std;
  39.  
  40. //----------------------------------------------------------------------
  41. PETProjDataMgr* PETProjDataMgr::m_Instance = 0;
  42.  
  43. //----------------------------------------------------------------------
  44. PETProjDataMgr* PETProjDataMgr::GetInstance()
  45. {
  46.   if( !m_Instance ) {
  47.     m_Instance = new PETProjDataMgr;
  48.   }
  49.  
  50.   return m_Instance;
  51.  
  52. }
  53.  
  54. //-----------------------------------------------------------------------
  55. PETProjDataMgr::PETProjDataMgr()
  56. {
  57.  
  58.   /*
  59. // ARGUMENTS:
  60. " cout << " -------------------------- \n"
  61. " Arguments convention: \n"
  62. "     -a Axial FOV (mm), <theDist_axial=100.0> \n"
  63. "     -d Diameter Transaxial FOV (mm), <m_RingDiameter=300.0> \n"
  64. "     -i Type of the input file (by default: 0 = Arce_binary), <typeINfile=0> \n"
  65. " \n"
  66. "     -m Maximum ring difference (by default: -1 = m_NOfPlanes), <m_MaxRingDifferenceiff> \n"
  67. "     -n Name of output file, <m_Filename> \n"
  68. "     -p Axial number of planes, <m_NOfPlanes> \n"
  69. "     -r Number of bins, \"distancias\", <m_NOfBins> \n"
  70.   // -s Span                                         TO DO:span !!!!!
  71. "     -t Number of angular views, \"direcciones\", <m_NOfAngles> \n"
  72. "     -v Verbosity (by default: 0=silent, 3=debug), <verbos> \n"
  73. "     -x Maximum number of coincidences to be stored (by default: -1 = no limit), <Max_Coinci> \n"
  74. "     -o Output type (by default: 0 = mcc Interfile, 1 = STIR Interfile), <OutType> \n"
  75.  
  76. " \n"
  77. " PET Reconstruction. CIEMAT 2009-11 \n"
  78. " mario.canadas@ciemat.es \n"
  79.  " -------------------------- \n";
  80.   */
  81.   m_rnd= new TRandom3();
  82.   m_AxialDistance = (9-1)*2.25; // Axial pixel dimension*NOfPlanes  
  83.   m_RingDiameter = 120.0; // notranji premer peta  
  84.   m_NOfPlanes = 9;        // stevilo ravnin
  85.   m_NOfBins = 128;        // stevilo binov v razdalji
  86.   m_nch =  128;            // stevilo padov okoli in okoli
  87.   m_NOfAngles = TMath::Abs(m_nch)/2;       // stevilo kotov = stevilo padov okoli in okoli /2
  88.  
  89.   m_MaxRingDifference = -1;
  90.   //m_MaxRingDifference = 3; // najvecja razdalja med padi
  91.   // toDo:  theSpan = int(GmParameterMgr::GetInstance()->GetNumericValue("PET:ProjData:Span",1));
  92.  
  93.  
  94.   m_OutFormat = 1; // 1.. projections 0 .. image
  95.   m_Filename = "sino3D";
  96.   m_Debug=1;
  97.  
  98.   if (m_MaxRingDifference==-1) m_MaxRingDifference=m_NOfPlanes-1;
  99.  
  100.  
  101.   m_TotalAxialPlanes=m_NOfPlanes*m_NOfPlanes;
  102.   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
  103.  
  104.   /*--- Initialize sino3D ---*/
  105.   m_projections = new SINO_TYPE**[m_NOfBins];
  106.   for(int i=0;i<m_NOfBins;i++){
  107.     m_projections[i] = new SINO_TYPE*[m_NOfAngles];
  108.     for(int j=0;j<m_NOfAngles;j++){
  109.       m_projections[i][j] = new SINO_TYPE[m_TotalAxialPlanes];  /// ! If m_OutFormat==1 (STIR output):Matrix size depends on the MAX_Ring_Difference
  110.       for(int k=0;k<m_TotalAxialPlanes;k++){
  111.         m_projections[i][j][k]=0;
  112.       }
  113.     }
  114.   }
  115.  
  116.   m_TotalProjectionCoincidences=0;
  117.   m_TotalCoincidences=0;
  118.   //OutputType = "pet";
  119. }
  120.  
  121. //-----------------------------------------------------------------------
  122. void PETProjDataMgr::SetProjection( int axialplane, TH2F * h)
  123. {
  124. for(int i=0;i<m_NOfBins;i++){
  125.   for(int j=0;j<m_NOfAngles;j++){
  126.     m_projections[i][j][axialplane]=h->GetBinContent(i+1,j+1);
  127.   }
  128. }
  129.  
  130. }
  131. //-----------------------------------------------------------------------
  132.  
  133. //-----------------------------------------------------------
  134. // from Gate
  135. //------------------------------------------------------------
  136. double ComputeSinogramS(double X1, double Y1,  double X2, double  Y2)
  137. {
  138.   double s;
  139.  
  140.   double denom = (Y1-Y2) * (Y1-Y2) + (X2-X1) * (X2-X1);
  141.  
  142.   if (denom!=0.) {
  143.     denom = sqrt(denom);
  144.     s = ( X1 * (Y2-Y1) + Y1 * (X1-X2)  ) / denom;
  145.   } else {
  146.     s = 0.;
  147.   }
  148.  
  149.   double theta;
  150.   if ((X1-X2)!=0.) {
  151.     theta=atan((X1-X2) /(Y1-Y2));
  152.   } else {
  153.     theta=3.1416/2.;
  154.   }
  155.   if ((theta > 0.) && ((X1-X2) > 0.)) s = -s;
  156.   if ((theta < 0.) && ((X1-X2) < 0.)) s = -s;
  157.   if ( theta < 0.) {
  158.     theta = theta+3.1416;
  159.     s = -s;
  160.   }
  161.   return s;
  162. }
  163.  
  164.  
  165. void PETProjDataMgr::AddEvent( const TVector3& pos1 , const TVector3& pos2)
  166. {
  167.   int z1_i, z2_i;
  168.     //for discretization on the crystal: int x1_i, x2_i, y1_i, y2_i;
  169.    
  170.   double z1_abs=pos1.z()+m_AxialDistance/2;
  171.   double z2_abs=pos2.z()+m_AxialDistance/2;
  172.   double a, b, phi, dis;
  173.   int phi_i, dis_i;
  174.   int ring_diff;
  175.  
  176.   double _PI=2*asin(1);
  177.  
  178.   m_TotalCoincidences++;
  179.  
  180.   z1_i=(int)(m_NOfPlanes* z1_abs/m_AxialDistance);  //round --> m_NOfPlanes+1 ...
  181.   z2_i=(int)(m_NOfPlanes* z2_abs/m_AxialDistance);
  182.  
  183.   // control; if z_i out of range: return
  184.  
  185.   if ( (pos1.x()==pos2.x()) && (pos1.y()==pos2.y()) ) {
  186. #ifndef GAMOS_NO_VERBOSE
  187.     if( m_Debug ) {
  188.       cout << "PETProjDataMgr::AddEvent:WARNING! Event_1 == Event_2  ; x= " << pos2.x() << ", y= " << pos2.y() << ", z= " << pos2.z() << endl;
  189.     }
  190. #endif
  191.     return;
  192.   }
  193.  
  194.   if ( (z1_i<0) || (z2_i<0) || (z1_i>= m_NOfPlanes) || (z2_i>= m_NOfPlanes) ) {
  195. #ifndef GAMOS_NO_VERBOSE
  196.     if( m_Debug ) {
  197.       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;
  198.     }
  199. #endif
  200.     return;
  201.   }
  202.  
  203.   ring_diff = (int)fabs(z1_i-z2_i);
  204.  
  205.   // max ring difference; control:
  206.   if (ring_diff > m_MaxRingDifference) {
  207. #ifndef GAMOS_NO_VERBOSE
  208.     if( m_Debug ) {
  209.       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;
  210.     }
  211. #endif
  212.     return;
  213.   }
  214.  
  215.   a=(double)(pos2.y()- pos1.y());
  216.   b=(double)(pos2.x()- pos1.x());
  217.  
  218.   if (a==0.0){
  219.     phi=_PI*0.5;
  220.   }
  221.   else{
  222.     phi=atan(b/a);
  223.   }
  224.  
  225.   if (phi<0) phi = phi +_PI;
  226.  
  227.   dis=pos1.x()*cos(phi) - pos1.y()*sin(phi);
  228.   //dis=ComputeSinogramS(pos1.x(), pos1.y(), pos2.x(), pos2.x());
  229.   // control; transaxial FOV
  230.   if ( fabs(dis) > m_RingDiameter*0.5 ) {
  231. #ifndef GAMOS_NO_VERBOSE
  232.     if( m_Debug ) {
  233.       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;
  234.     }
  235. #endif
  236.     return;
  237.   }
  238.  
  239.   dis = dis + m_RingDiameter*0.5;
  240.  
  241.   // discret values:
  242.   phi_i=(int)round( (double)(m_NOfAngles-1)*phi/_PI );
  243.   dis_i=(int)round( (double)(m_NOfBins-1)*dis/(double)m_RingDiameter );
  244.  
  245.   if ((phi_i>=m_NOfAngles) || (dis_i>=m_NOfBins))  return;  // only possible "=" because 'round' check it..
  246.  
  247.   // OLD: (SRRB included) sino3D[dis_i][phi_i][ (z1_i+z2_i)+ring_diff*(m_NOfPlanes-1) ]++;
  248.  
  249.   int Zpos;
  250.  
  251.   if (m_OutFormat==0) {
  252.     Zpos = (z1_i*m_NOfPlanes + z2_i);
  253.   }
  254.   else{
  255.    
  256.     if (z1_i>=z2_i) {  // SIN Max Ring_Diff: Zpos= ( ((m_NOfPlanes-ring_diff)*(m_NOfPlanes-1-ring_diff))/2 + z2_i );
  257.      
  258.       Zpos= ( ((2*m_NOfPlanes-1 - m_MaxRingDifference - ring_diff)*(m_MaxRingDifference - ring_diff))/2 + z2_i);
  259.      
  260.     }else{
  261.       Zpos= ( (m_TotalAxialPlanes) - ((2*m_NOfPlanes-1 - m_MaxRingDifference - ring_diff +1)*(m_MaxRingDifference - ring_diff +1))/2  + z1_i );
  262.      
  263.     }
  264.   }
  265.  
  266.   m_projections[dis_i][phi_i][ Zpos ]++;
  267.   m_TotalProjectionCoincidences++;
  268.  
  269. #ifndef GAMOS_NO_VERBOSE
  270.   if( m_Debug >1) {
  271.     cout << "PETProjDataMgr::AddEvent: x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl;
  272.     cout << "PETProjDataMgr::AddEvent: Sinogram pos.: distance(s)= " << dis_i << ", angular view(phi)= " << phi_i << " ; Zpos=" << Zpos <<"; Segment (Ring diff.) = " << ring_diff << endl;
  273.   }
  274. #endif
  275.  
  276.  
  277. }
  278.  
  279. //-----------------------------------------------------------------------
  280. PETProjDataMgr::~PETProjDataMgr()
  281. {
  282.  int i,j;
  283.                
  284.         for(i=0;i<m_NOfBins;i++){
  285.                 for(j=0;j<m_NOfAngles;j++){
  286.                         free(m_projections[i][j]);
  287.                
  288.                 }
  289.                 free(m_projections[i]);
  290.        
  291.         }
  292.         free(m_projections);
  293.  
  294. }
  295.  
  296. /*  TO DO: call lm_to_sino3D program
  297. //-----------------------------------------------------------------------
  298. void PETIOMgr::ReadFile()
  299. {
  300.   if( !theFileIn ) OpenFileIn();
  301.  
  302.   PETOutput po;
  303.   G4bool bEof;
  304.   for(;;) {
  305.     po = ReadEvent( bEof );
  306.     // theFileIn->read(reinterpret_cast<char *>(&po),sizeof(PetOutput));
  307.     if(bDumpCout) PrintEvent(" PETOutput: ", po, bCartesian);
  308.     if( bEof ) break;
  309.   }
  310. }
  311. //-----------------------------------------------------------------------
  312. PETOutput PETProjDataMgr::ReadEvent( G4bool& bEof )
  313. {
  314.   if( theFileIn == 0 ){
  315.     G4Exception("PETIOMgr::ReadEvent, file not opened, call OpenFileIn() first ");
  316.   }
  317.  
  318.   PETOutput po;
  319.   fread (&po, sizeof(struct PETOutput),1,theFileIn);
  320.   if ( feof( theFileIn ) ) {
  321.     bEof = TRUE;
  322.   } else {
  323.     bEof = FALSE;
  324.   }
  325.  
  326.   return po;
  327.  
  328. }
  329. */
  330.  
  331.  
  332. //-----------------------------------------------------------------------
  333. void PETProjDataMgr::WriteInterfile()
  334. {
  335.  
  336.   char name_hv[512];
  337.   char name_v[512];
  338.  
  339.   if (m_OutFormat==0){
  340.    
  341.     strcpy(name_hv, m_Filename);
  342.     strcpy(name_v,m_Filename);
  343.     strcat(name_hv, ".hv");
  344.     strcat(name_v, ".v");
  345.    
  346.     fp=fopen(name_hv, "w");
  347.    
  348.     fprintf (fp, "!INTERFILE  := \n");
  349.     fprintf (fp, "name of data file := %s\n", name_v);
  350.     fprintf (fp, "!GENERAL DATA := \n");
  351.     fprintf (fp, "!GENERAL IMAGE DATA :=\n");
  352.     fprintf (fp, "!type of data := tomographic\n");
  353.     fprintf (fp, "!version of keys := 3.3\n");
  354.     fprintf (fp, "!data offset in bytes := 0\n");
  355.     fprintf (fp, "imagedata byte order := littleendian\n");
  356.     fprintf (fp, "!PET STUDY (General) :=\n");
  357.     fprintf (fp, "!PET data type := 3D-Sinogram\n");
  358.     fprintf (fp, "process status := Reconstructed\n");
  359.     fprintf (fp, "!number format := unsigned short\n");
  360.     fprintf (fp, "!number of bytes per pixel := 2\n");
  361.     fprintf (fp, "number of dimensions := 3\n");
  362.     fprintf (fp, "matrix axis label [1] := x\n");
  363.     fprintf (fp, "!matrix size [1] := %i\n",m_NOfBins);
  364.     fprintf (fp, "scaling factor (mm/pixel) [1] := %f\n",(float)(m_RingDiameter/(m_NOfBins-1.0)));
  365.    
  366.     fprintf (fp, "matrix axis label [2] := y\n");
  367.     fprintf (fp, "!matrix size [2] := %i\n",m_NOfAngles);
  368.    
  369.     fprintf (fp, "scaling factor (degree/pixel) [2] := %f\n",(float)(360./m_NOfAngles));
  370.    
  371.     fprintf (fp, "matrix axis label [3] := z\n");
  372.     fprintf (fp, "!matrix size [3] := %i\n",m_NOfPlanes*m_NOfPlanes);
  373.     fprintf (fp, "scaling factor (mm/pixel) [3] := %f\n",(float)(m_AxialDistance/(m_NOfPlanes-1.0)));
  374.    
  375.     fprintf (fp, "number of slices := %i\n",m_NOfPlanes*m_NOfPlanes);
  376.     fprintf (fp, "number of time frames := 1\n");
  377.     fprintf (fp, "image scaling factor[1] := 1\n");
  378.     fprintf (fp, "data offset in bytes[1] := 0\n");
  379.     fprintf (fp, "quantification units := 1\n");
  380.     fprintf (fp, "!END OF INTERFILE := \n");
  381.    
  382.     fclose(fp);
  383.     //(size_t)(m_NOfBins*m_NOfAngles*m_NOfPlanes*m_NOfPlanes);
  384.    
  385.   }else{
  386.    
  387.     strcpy(name_hv, m_Filename);
  388.     strcpy(name_v,m_Filename);
  389.    
  390.     strcat(name_hv, ".hs"); // STIR extension: .hs .s
  391.     strcat(name_v, ".s");
  392.    
  393.     fp=fopen(name_hv, "w");
  394.    
  395.     fprintf (fp, "!INTERFILE  := \n");
  396.     fprintf (fp, "name of data file := %s\n",name_v);
  397.     fprintf (fp, "!GENERAL DATA := \n");
  398.     fprintf (fp, "!GENERAL IMAGE DATA :=\n");
  399.     fprintf (fp, "!type of data := PET\n");
  400.     //    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
  401.     //    fprintf (fp, "!data offset in bytes := 0\n");
  402.     fprintf (fp, "imagedata byte order := littleendian\n");
  403.     fprintf (fp, "!PET STUDY (General) :=\n");
  404.     fprintf (fp, "!PET data type := Emission\n");
  405.     fprintf (fp, "applied corrections := {arc correction}\n");        // {none}\n");
  406.     //    fprintf (fp, "process status := Reconstructed\n");
  407.     fprintf (fp, "!number format := unsigned integer\n");
  408.     fprintf (fp, "!number of bytes per pixel := 2\n");
  409.  
  410.     fprintf (fp, "number of dimensions := 4\n");
  411.     fprintf (fp, "matrix axis label [4] := segment\n");
  412.     fprintf (fp, "!matrix size [4] := %i\n",m_MaxRingDifference*2 + 1);
  413.     //    fprintf (fp, "scaling factor (mm/pixel) [1] := %f\n",(float)(d_FOV/(m_NOfBins-1)));
  414.     fprintf (fp, "matrix axis label [3] := axial coordinate\n");
  415.     fprintf (fp, "!matrix size [3] := { ");
  416.     if (m_MaxRingDifference==0)
  417.       {
  418.         fprintf (fp, "%i}\n", m_NOfPlanes);
  419.       }else{
  420.         for(int m=m_NOfPlanes-m_MaxRingDifference;m<=m_NOfPlanes;m++)  fprintf (fp, "%i,", m);
  421.         for(int m=m_NOfPlanes-1;m>m_NOfPlanes-m_MaxRingDifference;m--)  fprintf (fp, "%i,", m);
  422.         fprintf (fp, "%i}\n", m_NOfPlanes-m_MaxRingDifference);
  423.       }
  424.     fprintf (fp, "matrix axis label [2] := view\n");
  425.     fprintf (fp, "!matrix size [2] := %i\n",m_NOfAngles);
  426.     fprintf (fp, "matrix axis label [1] := tangential coordinate\n");
  427.     fprintf (fp, "!matrix size [1] := %i\n",m_NOfBins);  
  428.  
  429.     fprintf (fp, "minimum ring difference per segment := {");     // TO DO :  add SPAN (m_MaxRingDifferenceiff per seg. variable)
  430.     fprintf (fp, "%i", -m_MaxRingDifference);
  431.     for(int m=-m_MaxRingDifference+1;m<=m_MaxRingDifference;m++)  fprintf (fp, ",%i", m);
  432.     fprintf (fp, "}\n");
  433.     fprintf (fp, "maximum ring difference per segment := {");     // TO DO :  add SPAN (m_MaxRingDifferenceiff per seg. variable)
  434.     fprintf (fp, "%i", -m_MaxRingDifference);
  435.     for(int m=-m_MaxRingDifference+1;m<=m_MaxRingDifference;m++)  fprintf (fp, ",%i", m);
  436.     fprintf (fp, "}\n");
  437.  
  438.     fprintf (fp, "inner ring diameter (cm) := %f\n", m_RingDiameter/10);   // STIR Required parameter, now assigned to FOV (not detectors)
  439.     fprintf (fp, "average depth of interaction (cm) := 0.0001\n");
  440.     fprintf (fp, "default bin size (cm) := %f\n",0.1*((float)m_RingDiameter/((float)m_NOfBins-1.0)) );
  441.     fprintf (fp, "number of rings := %i\n",m_NOfPlanes );    
  442.     fprintf (fp, "distance between rings (cm) := %f\n", 0.1*((float)m_AxialDistance/(float)(m_NOfPlanes-1)) );  // Axial pixel dimension
  443.    
  444.     fprintf (fp, "number of detectors per ring := %i\n",m_NOfAngles*2 );  
  445.    //    fprintf (fp, "number of slices := %i\n",m_NOfPlanes*m_NOfPlanes);
  446.     fprintf (fp, "number of time frames := 1\n");
  447.     fprintf (fp, "image scaling factor[1] := 1\n");
  448.     fprintf (fp, "data offset in bytes[1] := 0\n");
  449.     fprintf (fp, "quantification units := 1\n");
  450.     fprintf (fp, "!END OF INTERFILE := \n");
  451.    
  452.     fclose(fp);
  453.  
  454.   }
  455.   m_Buffer = (SINO_TYPE*) malloc( m_NOfBins*m_NOfAngles*m_TotalAxialPlanes*sizeof(SINO_TYPE));
  456.  
  457.   long unsigned int cont=0;
  458.   int i,j,k;
  459.  
  460.   for(k=0;k<m_TotalAxialPlanes;k++){
  461.     for(j=0;j<m_NOfAngles;j++){
  462.       for(i=0;i<m_NOfBins;i++){
  463.         m_Buffer[cont]=m_projections[i][j][k];
  464.         cont++;
  465.       }
  466.     }
  467.   }
  468.  
  469.   fp=fopen(name_v, "wb");
  470.  
  471.   //cout << 4096*sizeof(SINO_TYPE) << endl;
  472.   int nb=fwrite(m_Buffer,1,m_NOfBins*m_NOfAngles*m_TotalAxialPlanes*sizeof(SINO_TYPE), fp);
  473.   fclose(fp);
  474.  
  475. #ifndef GAMOS_NO_VERBOSE
  476.   cout << "PETProjDataMgr::WriteInterfile: File name: " << m_Filename << endl;
  477.   cout << "PETProjDataMgr::WriteInterfile: Numer of bytes written: " << nb << endl;
  478.   cout << "PETProjDataMgr::WriteInterfile: Planes = " << m_NOfPlanes << "; bins = " << m_NOfBins << "; ang_views = " << m_NOfAngles << endl;
  479.   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);
  480.   cout << "... " << endl;
  481.  
  482.  cout << "PETProjDataMgr::WriteInterfile: Total  Coinci: " << m_TotalCoincidences << endl;
  483.  cout << "PETProjDataMgr::WriteInterfile: Sino3D Coinci: " << m_TotalProjectionCoincidences << endl;  
  484. #endif
  485.  
  486. }
  487.  
  488.  
  489. TVector3 PETProjDataMgr::Hits2Digits(const TVector3 &r){
  490.   if (!m_nch) return r;
  491.   float smear=0.5;
  492.  
  493.   if (m_nch<0) smear=m_rnd->Rndm();
  494.  
  495.   double angle = TMath::ATan2(r.X(),r.Y());  // vrne kot med -pi in pi
  496.   if (angle<0) angle+=TMath::TwoPi();
  497.  
  498.   angle= (int(angle/TMath::TwoPi()*TMath::Abs(m_nch))+smear)*TMath::TwoPi()/TMath::Abs(m_nch);
  499.   //(m_rnd->Rndm()-0.5)*m_AxialDistance;
  500.   return TVector3(sin(angle), cos(angle),0); // z coordinata ni cisto v redu
  501.  
  502. }
  503.  
  504. int PETProjDataMgr::FwdProject(double x,double y, double z, int nmax, TH1*h){
  505. TVector3 r(x,y,z);
  506. int h2d=h->InheritsFrom("TH2F");
  507. double tfac=m_RingDiameter*m_RingDiameter/4-r.Mag2();
  508. double rfac= m_AxialDistance/m_RingDiameter;
  509. for (int i=0;i<nmax;i++){
  510.  
  511.   double phi= m_rnd->Rndm()*TMath::Pi();  
  512.   TVector3 s(1,0,0);
  513.   s.SetPhi(phi);
  514.   double sign = (m_rnd->Rndm()>0.5)? 1 : 0;
  515.   double theta = TMath::ACos(m_rnd->Rndm()*rfac);
  516.   theta+=sign*TMath::Pi();
  517.  
  518.   s.SetTheta(theta);
  519.   double t=r*s;
  520.   TVector3 rx=r-t*s;
  521.  
  522.   double d=TMath::Sqrt(t*t+tfac);
  523.  
  524.   TVector3 r1=rx+d*s;
  525.   TVector3 r2=rx-d*s;
  526.  
  527.   //r1=Hits2Digits(r1);
  528.   //r2=Hits2Digits(r2);
  529.  
  530.   AddEvent( r1 , r2);
  531.   if (h!=NULL){
  532.     TVector3 s1=r2-r1;
  533.     double s1len= s1.Mag();
  534.     int niter=int (100*s1len/m_RingDiameter);
  535.     for (int j=0;j<niter;j++){
  536.       r2=r1+m_rnd->Rndm()*s1;
  537.       if (h2d) ((TH2F *) h)->Fill(r2.X(),r2.Y());
  538.       else ((TH3F *) h)->Fill(r2.X(),r2.Y(),r2.Z());
  539.     }
  540.   }
  541. }
  542. return 0;
  543. }
  544.  
  545. int PETProjDataMgr::FwdProject(TH2F *img, TH2F *h){
  546.  
  547. for (int i=0;i<img->GetNbinsX();i++) {
  548.    double x_=img->GetXaxis()->GetBinCenter( i+1 );
  549.    for (int j=0;j<img->GetNbinsY();j++) {
  550.        double y_=img->GetYaxis()->GetBinCenter( j+1 );
  551.        double density= img->GetBinContent(i+1,j+1);
  552.        if (density>0) FwdProject(x_,y_,m_AxialDistance*(m_rnd->Rndm()-0.5), density,h);
  553.    }
  554. }
  555. return 0;
  556. }
  557.  
  558. int PETProjDataMgr::FwdProject(TH3F *img, TH3F *h){
  559.  
  560. for (int i=0;i<img->GetNbinsX();i++) {
  561.    double x_=img->GetXaxis()->GetBinCenter( i+1 );
  562.    for (int j=0;j<img->GetNbinsY();j++) {
  563.       double y_=img->GetYaxis()->GetBinCenter( j+1 );
  564.       for (int k=0;k<img->GetNbinsZ();k++) {
  565.          double z_=img->GetZaxis()->GetBinCenter( k+1 );
  566.          double density= img->GetBinContent(i+1,j+1,k+1);
  567.          if (density>0) FwdProject(x_,y_,z_, density,h);
  568.       }
  569.    }
  570. }
  571. return 0;
  572. }
  573.  
  574. TH2F *PETProjDataMgr::Phantom(int kaj){
  575. TH2F *img= new TH2F("img","Original Image",100,-50,50,100,-50,50);
  576.  
  577. // izberi sliko 0: kroglice, 1: point source 2: central ball
  578. switch (kaj){
  579.  
  580. case 0:
  581. for (int i=0;i<img->GetNbinsX();i++) {
  582.    for (int j=0;j<img->GetNbinsY();j++) {
  583.        double x_=img->GetXaxis()->GetBinCenter( i+1 );
  584.        double y_=img->GetYaxis()->GetBinCenter( j+1 );
  585.        double density=1000;
  586.        if ((x_*x_+y_*y_)<6) img->SetBinContent(i+1,j+1,density);
  587.        
  588.        density=500; if ((x_-25)*(x_-25)+y_*y_<12) img->SetBinContent(i+1,j+1,density);
  589.        density=2000; if ((y_-25)*(y_-25)+x_*x_<2) img->SetBinContent(i+1,j+1,density);
  590.    }
  591. }
  592. break;
  593.  
  594. case 2:
  595. for (int i=0;i<img->GetNbinsX();i++) {
  596.    for (int j=0;j<img->GetNbinsY();j++) {
  597.        double x_=img->GetXaxis()->GetBinCenter( i+1 );
  598.        double y_=img->GetYaxis()->GetBinCenter( j+1 );
  599.        double density=1000;
  600.        if ((x_*x_+y_*y_)<12.5) img->SetBinContent(i+1,j+1,density);
  601.    }
  602. }
  603. break;
  604.  
  605. case 1:
  606. img->Fill(25,25,10000);
  607. break;
  608.  
  609. }
  610.  
  611. return img;
  612. }
  613.