Subversion Repositories f9daq

Rev

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