WebSVN - f9daq - Rev 264 - /cvi/apps/petdemo/PETProjDataMgr.c

//
// ../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"

static struct PETProjDataMgr mgr;

void SetNPlanes(int n) {
mgr.m_NOfPlanes=n;
};
void SetNBin(int n) {
mgr.m_NOfBins=n;
};
void SetNAng(int n) {
mgr.m_NOfAngles=n;
};
void SetDebug(int n) {
mgr.m_Debug=n;
};
void SetNumberOfChannels(int n) {
mgr.m_nch=n;
};

void SetRingDiameter( double x) {
mgr.m_RingDiameter = x;
};
void SetAxialDistance( double x) {
mgr.m_AxialDistance = x;
};

//----------------------------------------------------------------------

//----------------------------------------------------------------------
struct PETProjDataMgr *GetInstance() {

return &mgr;

}

//-----------------------------------------------------------------------
int PETProjDataMgrInit() {

/*
// 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";
*/

mgr.m_AxialDistance = (9-1)*2.25; // Axial pixel dimension*NOfPlanes
mgr.m_RingDiameter = 100.0; // notranji premer peta
mgr.m_NOfPlanes = 9; // stevilo ravnin
mgr.m_NOfBins = 128; // stevilo binov v razdalji
mgr.m_nch = 128; // stevilo padov okoli in okoli
mgr.m_NOfAngles = fabs(mgr.m_nch)/2; // stevilo kotov = stevilo padov okoli in okoli /2

mgr.m_MaxRingDifference = -1;
//mgr.m_MaxRingDifference = 3; // najvecja razdalja med padi
// toDo: theSpan = int(GmParameterMgr::GetInstance()->GetNumericValue("PET:ProjData:Span",1));

mgr.m_OutFormat = 1; // 1.. projections 0 .. image

mgr.m_Debug=1;

if (mgr.m_MaxRingDifference==-1) mgr.m_MaxRingDifference=mgr.m_NOfPlanes-1;

mgr.m_TotalAxialPlanes=mgr.m_NOfPlanes*mgr.m_NOfPlanes;
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

/*--- Initialize sino3D ---*/
mgr.m_projections = (SINO_TYPE ** *) malloc(mgr.m_NOfBins*sizeof(SINO_TYPE **));
for(int i=0; i<mgr.m_NOfBins; i++) {
mgr.m_projections[i] = (SINO_TYPE **) malloc(mgr.m_NOfAngles*sizeof(SINO_TYPE *));
for(int j=0; j<mgr.m_NOfAngles; j++) {
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
for(int k=0; k<mgr.m_TotalAxialPlanes; k++) {
mgr.m_projections[i][j][k]=0;
}
}
}

mgr.m_TotalProjectionCoincidences=0;
mgr.m_TotalCoincidences=0;
//OutputType = "pet";
return 0;
}

//-----------------------------------------------------------------------
void SetProjection( int axialplane, int id) {
for(int i=0; i<mgr.m_NOfBins; i++) {
for(int j=0; j<mgr.m_NOfAngles; j++) {
mgr.m_projections[i][j][axialplane]=H2D_GetBinContent(id,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 AddEvent( const HVector3 pos1 , const HVector3 pos2) {
int z1_i, z2_i;
//for discretization on the crystal: int x1_i, x2_i, y1_i, y2_i;

double z1_abs=pos1.x[2]+mgr.m_AxialDistance/2;
double z2_abs=pos2.x[2]+mgr.m_AxialDistance/2;
double a, b, phi, dis;
int phi_i, dis_i;
int ring_diff;

double _PI=2*asin(1);

mgr.m_TotalCoincidences++;

z1_i=(int)(mgr.m_NOfPlanes* z1_abs/mgr.m_AxialDistance); //round --> mgr.m_NOfPlanes+1 ...
z2_i=(int)(mgr.m_NOfPlanes* z2_abs/mgr.m_AxialDistance);

// control; if z_i out of range: return

if ( (pos1.x[0]==pos2.x[0]) && (pos1.x[1]==pos2.x[1]) ) {
#ifndef GAMOS_NO_VERBOSE
if( mgr.m_Debug ) {
printf( "AddEvent:WARNING! Event_1 == Event_2 ; x= %f y= %f z= %f\n", pos2.x[0], pos2.x[1], pos2.x[2] );
}
#endif
return;
}

if ( (z1_i<0) || (z2_i<0) || (z1_i>= mgr.m_NOfPlanes) || (z2_i>= mgr.m_NOfPlanes) ) {
#ifndef GAMOS_NO_VERBOSE
if( mgr.m_Debug ) {
printf("PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Axial):");
printf( "x1= %f y1= %f z1= %f ", pos1.x[0], pos1.x[1], pos1.x[2] );
printf( "x2= %f y2= %f z2= %f \n", pos2.x[0], pos2.x[1], pos2.x[2] );
}
#endif
return;
}

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

// max ring difference; control:
if (ring_diff > mgr.m_MaxRingDifference) {
#ifndef GAMOS_NO_VERBOSE
if( mgr.m_Debug ) {
printf("PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Max. Ring Diff.): %f>%f",ring_diff , mgr.m_MaxRingDifference );
printf( "x1= %f y1= %f z1= %f ", pos1.x[0], pos1.x[1], pos1.x[2] );
printf( "x2= %f y2= %f z2= %f \n", pos2.x[0], pos2.x[1], pos2.x[2] );
}
#endif
return;
}

a=(double)(pos2.x[1]- pos1.x[1]);
b=(double)(pos2.x[0]- pos1.x[0]);

if (a==0.0) {
phi=_PI*0.5;
} else {
phi=atan(b/a);
}

if (phi<0) phi = phi +_PI;

dis=pos1.x[0]*cos(phi) - pos1.x[1]*sin(phi);
//dis=ComputeSinogramS(pos1.x(), pos1.y(), pos2.x(), pos2.x());
// control; transaxial FOV
if ( fabs(dis) > mgr.m_RingDiameter*0.5 ) {
#ifndef GAMOS_NO_VERBOSE
if( mgr.m_Debug ) {
printf("PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Transaxial):" );
printf( "x1= %f y1= %f z1= %f ", pos1.x[0], pos1.x[1], pos1.x[2] );
printf( "x2= %f y2= %f z2= %f \n", pos2.x[0], pos2.x[1], pos2.x[2] );
}
#endif
return;
}

dis = dis + mgr.m_RingDiameter*0.5;

// discret values:
phi_i=RoundRealToNearestInteger( (double)(mgr.m_NOfAngles-1)*phi/_PI );
dis_i=RoundRealToNearestInteger( (double)(mgr.m_NOfBins-1)*dis/(double)mgr.m_RingDiameter );

if ((phi_i>=mgr.m_NOfAngles) || (dis_i>=mgr.m_NOfBins)) return; // only possible "=" because 'round' check it..

// OLD: (SRRB included) sino3D[dis_i][phi_i][ (z1_i+z2_i)+ring_diff*(mgr.m_NOfPlanes-1) ]++;

int Zpos;

if (mgr.m_OutFormat==0) {
Zpos = (z1_i*mgr.m_NOfPlanes + z2_i);
} else {

if (z1_i>=z2_i) { // SIN Max Ring_Diff: Zpos= ( ((mgr.m_NOfPlanes-ring_diff)*(mgr.m_NOfPlanes-1-ring_diff))/2 + z2_i );

Zpos= ( ((2*mgr.m_NOfPlanes-1 - mgr.m_MaxRingDifference - ring_diff)*(mgr.m_MaxRingDifference - ring_diff))/2 + z2_i);

} else {
Zpos= ( (mgr.m_TotalAxialPlanes) - ((2*mgr.m_NOfPlanes-1 - mgr.m_MaxRingDifference - ring_diff +1)*(mgr.m_MaxRingDifference - ring_diff +1))/2 + z1_i );

}
}

mgr.m_projections[dis_i][phi_i][ Zpos ]++;
mgr.m_TotalProjectionCoincidences++;

#ifndef GAMOS_NO_VERBOSE
if( mgr.m_Debug >1) {
printf("PETProjDataMgr::AddEvent:");
printf( "x1= %f y1= %f z1= %f ", pos1.x[0], pos1.x[1], pos1.x[2] );
printf( "x2= %f y2= %f z2= %f \n", pos2.x[0], pos2.x[1], pos2.x[2] );
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 );
}
#endif

}

//-----------------------------------------------------------------------
void PETProjDataMgrFree() {
int i,j;
if (mgr.m_projections){
for(i=0; i<mgr.m_NOfBins; i++) {
for(j=0; j<mgr.m_NOfAngles; j++) {
if (mgr.m_projections[i][j]) free(mgr.m_projections[i][j]);
mgr.m_projections[i][j]= NULL;
}
if (mgr.m_projections[i]) free(mgr.m_projections[i]);
mgr.m_projections[i]= NULL;

}
free(mgr.m_projections);
mgr.m_projections= NULL;
}

}

//-----------------------------------------------------------------------
void WriteInterfile(const char * Filename ) {

char name_hv[512];
char name_v[512];
FILE *fp;
if (mgr.m_OutFormat==0) {

strcpy(name_hv, Filename);
strcpy(name_v, 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",mgr.m_NOfBins);
fprintf (fp, "scaling factor (mm/pixel) [1] := %f\n",(float)(mgr.m_RingDiameter/(mgr.m_NOfBins-1.0)));

fprintf (fp, "matrix axis label [2] := y\n");
fprintf (fp, "!matrix size [2] := %i\n",mgr.m_NOfAngles);

fprintf (fp, "scaling factor (degree/pixel) [2] := %f\n",(float)(360./mgr.m_NOfAngles));

fprintf (fp, "matrix axis label [3] := z\n");
fprintf (fp, "!matrix size [3] := %i\n",mgr.m_NOfPlanes*mgr.m_NOfPlanes);
fprintf (fp, "scaling factor (mm/pixel) [3] := %f\n",(float)(mgr.m_AxialDistance/(mgr.m_NOfPlanes-1.0)));

fprintf (fp, "number of slices := %i\n",mgr.m_NOfPlanes*mgr.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)(mgr.m_NOfBins*mgr.m_NOfAngles*mgr.m_NOfPlanes*mgr.m_NOfPlanes);

} else {

strcpy(name_hv, Filename);
strcpy(name_v, 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",mgr.m_MaxRingDifference*2 + 1);
// fprintf (fp, "scaling factor (mm/pixel) [1] := %f\n",(float)(d_FOV/(mgr.m_NOfBins-1)));
fprintf (fp, "matrix axis label [3] := axial coordinate\n");
fprintf (fp, "!matrix size [3] := { ");
if (mgr.m_MaxRingDifference==0) {
fprintf (fp, "%i}\n", mgr.m_NOfPlanes);
} else {
int m;
for(m=mgr.m_NOfPlanes-mgr.m_MaxRingDifference; m<=mgr.m_NOfPlanes; m++) fprintf (fp, "%i,", m);
for(m=mgr.m_NOfPlanes-1; m>mgr.m_NOfPlanes-mgr.m_MaxRingDifference; m--) fprintf (fp, "%i,", m);
fprintf (fp, "%i}\n", mgr.m_NOfPlanes-mgr.m_MaxRingDifference);
}
fprintf (fp, "matrix axis label [2] := view\n");
fprintf (fp, "!matrix size [2] := %i\n",mgr.m_NOfAngles);
fprintf (fp, "matrix axis label [1] := tangential coordinate\n");
fprintf (fp, "!matrix size [1] := %i\n",mgr.m_NOfBins);

fprintf (fp, "minimum ring difference per segment := {"); // TO DO : add SPAN (mgr.m_MaxRingDifferenceiff per seg. variable)
fprintf (fp, "%i", -mgr.m_MaxRingDifference);
int m;
for(m=-mgr.m_MaxRingDifference+1; m<=mgr.m_MaxRingDifference; m++) fprintf (fp, ",%i", m);
fprintf (fp, "}\n");
fprintf (fp, "maximum ring difference per segment := {"); // TO DO : add SPAN (mgr.m_MaxRingDifferenceiff per seg. variable)
fprintf (fp, "%i", -mgr.m_MaxRingDifference);
for(m=-mgr.m_MaxRingDifference+1; m<=mgr.m_MaxRingDifference; m++) fprintf (fp, ",%i", m);
fprintf (fp, "}\n");

fprintf (fp, "inner ring diameter (cm) := %f\n", mgr.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)mgr.m_RingDiameter/((float)mgr.m_NOfBins-1.0)) );
fprintf (fp, "number of rings := %i\n",mgr.m_NOfPlanes );
fprintf (fp, "distance between rings (cm) := %f\n", 0.1*((float)mgr.m_AxialDistance/(float)(mgr.m_NOfPlanes-1)) ); // Axial pixel dimension

fprintf (fp, "number of detectors per ring := %i\n",mgr.m_NOfAngles*2 );
// fprintf (fp, "number of slices := %i\n",mgr.m_NOfPlanes*mgr.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);

}
int mbsize = mgr.m_NOfBins*mgr.m_NOfAngles*mgr.m_TotalAxialPlanes*sizeof(SINO_TYPE);
mgr.m_Buffer = (SINO_TYPE *) malloc( mbsize );

long unsigned int cont=0;
int i,j,k;

for(k=0; k<mgr.m_TotalAxialPlanes; k++) {
for(j=0; j<mgr.m_NOfAngles; j++) {
for(i=0; i<mgr.m_NOfBins; i++) {
mgr.m_Buffer[cont]=mgr.m_projections[i][j][k];
cont++;
}
}
}

fp=fopen(name_v, "wb");

//printf(4096*sizeof(SINO_TYPE) );
int nb=fwrite(mgr.m_Buffer,1,mbsize, fp);
fclose(fp);

#ifndef GAMOS_NO_VERBOSE
printf("PETProjDataMgr::WriteInterfile: File name: %s\n", Filename );
printf("PETProjDataMgr::WriteInterfile: Numer of bytes written: %d\n" , nb );
printf("PETProjDataMgr::WriteInterfile: Planes = %d bins = %d ang_views = %d\n", mgr.m_NOfPlanes, mgr.m_NOfBins, mgr.m_NOfAngles );
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);
printf("... " );

printf("PETProjDataMgr::WriteInterfile: Total Coinci: %d\n" , mgr.m_TotalCoincidences );
printf("PETProjDataMgr::WriteInterfile: Sino3D Coinci: %d\n" ,mgr.m_TotalProjectionCoincidences );
#endif

}

double Mag2(HVector3 x) {
return x.x[0]*x.x[0]+x.x[1]*x.x[1]+x.x[2]*x.x[2];
}

double Mag(HVector3 x) {
return sqrt(Mag2(x));
}

void SetPhi(HVector3 *r,double phi) {

}
void SetTheta(HVector3 *r,double phi) {

}

double Product(HVector3 a, HVector3 b) {
double c=0;
for (int i=0; i<3; i++) c += a.x[i] * b.x[i];
return c;
}

HVector3 Multiply(double a, HVector3 b) {
HVector3 c;
for (int i=0; i<3; i++) c.x[i]= a* b.x[i];
return c;
}

HVector3 Add(HVector3 a, HVector3 b) {
HVector3 c;
for (int i=0; i<3; i++) c.x[i]= a.x[i] + b.x[i];
return c;
}

HVector3 Subtract(HVector3 a, HVector3 b) {
HVector3 c;
for (int i=0; i<3; i++) c.x[i]= a.x[i] - b.x[i];
return c;
}

HVector3 Hits2Digits(const HVector3 r) {
if (!mgr.m_nch) return r;
float smear=0.5;

if (mgr.m_nch<0) smear=Random(0,1);

double angle = atan2(r.x[0],r.x[1]); // vrne kot med -pi in pi
double twopi=2*Pi();
if (angle<0) angle+=twopi;

angle= ((int)(angle/twopi*fabs(mgr.m_nch))+smear)*twopi/fabs(mgr.m_nch);
//(mgr.m_rnd->Rndm()-0.5)*mgr.m_AxialDistance;
HVector3 x;
x.x[0]=sin(angle);
x.x[1]=cos(angle);
x.x[2]=0;
return x; // z coordinata ni cisto v redu

}

int FwdProject(double x,double y, double z, int nmax, int h) {
HVector3 r;
r.x[0]=x;
r.x[1]=y;
r.x[2]=z;
int h2d=1;
double tfac=mgr.m_RingDiameter*mgr.m_RingDiameter/4-Mag2(r);
double rfac= mgr.m_AxialDistance/mgr.m_RingDiameter;
for (int i=0; i<nmax; i++) {

double phi= Random(0,Pi());
HVector3 s;
s.x[0]=1;
s.x[1]=0;
s.x[2]=0;
SetPhi(&s,phi);
double sign = (Random(0,1)>0.5)? 1 : 0;
double theta = acos(Random(0,rfac));
theta+=sign*Pi();

SetTheta(&s,theta);
double t=Product(r,s);
HVector3 rx=Add(r,Multiply(-t,s));

double d=sqrt(t*t+tfac);

HVector3 r1= Add(rx,Multiply(d,s));
HVector3 r2= Add(rx,Multiply(-d,s));

//r1=Hits2Digits(r1);
//r2=Hits2Digits(r2);

HVector3 s1=Subtract(r2,r1);
double s1len= Mag(s1);
int niter=(int) (100*s1len/mgr.m_RingDiameter);
for (int j=0; j<niter; j++) {
r2=Add(r1,Multiply(Random(0,1),s1));
if (h2d) H2D_Fill(h,r2.x[0],r2.x[1],1);
else H3D_Fill(h,r2.x[0],r2.x[1],r2.x[2],1);
}
}
return 0;
}

int FwdProject2d(int img, int h) {

for (int i=0; i<H2D_GetNbinsX(img); i++) {
double x_=H2D_GetXBinCenter( img, i+1 );
for (int j=0; j<H2D_GetNbinsY(img); j++) {
double y_=H2D_GetYBinCenter(img, j+1 );
double density= H2D_GetBinContent(img,i+1,j+1);
if (density>0) FwdProject(x_,y_,mgr.m_AxialDistance*(Random(-0.5,0.5)), density,h);
}
}
return 0;
}

int FwdProject3d(int img, int h) {

for (int i=0; i<H3D_GetNbinsX(img); i++) {
double x_=H3D_GetXBinCenter( img, i+1 );
for (int j=0; j<H3D_GetNbinsY(img); j++) {
double y_=H3D_GetYBinCenter(img, j+1 );
for (int k=0; k<H3D_GetNbinsZ(img); k++) {
double z_=H3D_GetZBinCenter( img, k+1 );
double density= H3D_GetBinContent(img, i+1,j+1,k+1);
if (density>0) FwdProject(x_,y_,z_, density,h);
}
}
}
return 0;
}

int Phantom(int kaj) {
int img = H2D_Init(1, "img","Original Image",100,-50,50,100,-50,50);
int i=0;
// izberi sliko 0: kroglice, 1: point source 2: central ball
switch (kaj) {

case 0:
for (i=0; i<H2D_GetNbinsX(img); i++) {
for (int j=0; j<H2D_GetNbinsY(img); j++) {
double x_=H2D_GetXBinCenter( img, i+1 );
double y_=H2D_GetYBinCenter( img, j+1 );
double density=1000;
if ((x_*x_+y_*y_)<6) H2D_SetBinContent(img, i+1,j+1,density);

density=500;
if ((x_-25)*(x_-25)+y_*y_<12) H2D_SetBinContent(img,i+1,j+1,density);
density=2000;
if ((y_-25)*(y_-25)+x_*x_<2) H2D_SetBinContent(img,i+1,j+1,density);
}
}
break;

case 2:
for (i=0; i<H2D_GetNbinsX(img); i++) {
for (int j=0; j<H2D_GetNbinsY(img); j++) {
double x_=H2D_GetXBinCenter( img, i+1 );
double y_=H2D_GetYBinCenter( img, j+1 );
double density=1000;
if ((x_*x_+y_*y_)<12.5) H2D_SetBinContent(img, i+1,j+1,density);
}
}
break;

case 1:
H2D_Fill(img, 25,25,10000);
break;

}

return img;
}

Subversion Repositories f9daq

(root)/cvi/apps/petdemo/PETProjDataMgr.c - Rev 264