Subversion Repositories f9daq

Rev

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

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