Subversion Repositories f9daq

Rev

Blame | 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*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;
}