Subversion Repositories f9daq

//

// ../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*NOfPlanes

  m_RingDiameter = 120.0; // notranji premer peta

  m_NOfPlanes = 9;        // stevilo ravnin

  m_NOfBins = 128;        // stevilo binov v razdalji

  m_nch =  128;            // stevilo padov okoli in okoli

  m_NOfAngles = TMath::Abs(m_nch)/2;       // stevilo kotov = stevilo padov okoli in okoli /2

  m_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 .. image

  m_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_Difference

      for(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: return

  if ( (pos1.x()==pos2.x()) && (pos1.y()==pos2.y()) ) {

#ifndef GAMOS_NO_VERBOSE

    if( m_Debug ) {

      cout << "PETProjDataMgr::AddEvent:WARNING! Event_1 == Event_2  ; x= " << pos2.x() << ", y= " << pos2.y() << ", z= " << pos2.z() << endl;

    }

#endif

    return;

  }

  if ( (z1_i<0) || (z2_i<0) || (z1_i>= m_NOfPlanes) || (z2_i>= m_NOfPlanes) ) {

#ifndef GAMOS_NO_VERBOSE

    if( 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;

    }

#endif

    return;

  }

  ring_diff = (int)fabs(z1_i-z2_i);

  // max ring difference; control:

  if (ring_diff > m_MaxRingDifference) {

#ifndef GAMOS_NO_VERBOSE

    if( 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;

    }

#endif

    return;

  }

  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 FOV

  if ( fabs(dis) > m_RingDiameter*0.5 ) {

#ifndef GAMOS_NO_VERBOSE

    if( 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;

    }

#endif

    return;

  }

  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_VERBOSE

  if( 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 .s

    strcat(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 dimension

    fprintf (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_VERBOSE

  cout << "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 pi

  if (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 ball

switch (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;

}