CASToR  1.1
Tomographic Reconstruction (PET/SPECT)
 All Classes Files Functions Variables Typedefs Macros Groups Pages
iProjectorJoseph.cc
Go to the documentation of this file.
1 
2 /*
3  Implementation of class iProjectorJoseph
4 
5  - separators: done
6  - doxygen: done
7  - default initialization: done
8  - CASTOR_DEBUG:
9  - CASTOR_VERBOSE:
10 */
11 
18 #include "iProjectorJoseph.hh"
19 #include "sOutputManager.hh"
20 
21 #include <cmath>
22 
23 // =====================================================================
24 // ---------------------------------------------------------------------
25 // ---------------------------------------------------------------------
26 // =====================================================================
27 
28 iProjectorJoseph::iProjectorJoseph() : vProjector()
29 {
30  // This projector is not compatible with SPECT attenuation correction because
31  // the voxels contributing to the line are not strictly ordered with respect to
32  // their distance to point2 (due to interpolations at each plane crossed)
33  m_compatibleWithSPECTAttenuationCorrection = false;
34 
35  // Default pointers and parameters
36  mp_boundariesX = nullptr;
37  mp_boundariesY = nullptr;
38  mp_boundariesZ = nullptr;
39  mp_maskPad = nullptr;
40  m_toleranceX = 0.;
41  m_toleranceY = 0.;
42  m_toleranceZ = 0.;
43 }
44 
45 // =====================================================================
46 // ---------------------------------------------------------------------
47 // ---------------------------------------------------------------------
48 // =====================================================================
49 
50 iProjectorJoseph::~iProjectorJoseph()
51 {
52  if ( mp_boundariesX )
53  {
54  delete[] mp_boundariesX;
55  mp_boundariesX = nullptr;
56  }
57 
58  if ( mp_boundariesY )
59  {
60  delete[] mp_boundariesY;
61  mp_boundariesY = nullptr;
62  }
63 
64  if ( mp_boundariesZ )
65  {
66  delete[] mp_boundariesZ;
67  mp_boundariesZ = nullptr;
68  }
69 
70  if ( mp_maskPad )
71  {
72  delete[] mp_maskPad;
73  mp_maskPad = nullptr;
74  }
75 }
76 
77 // =====================================================================
78 // ---------------------------------------------------------------------
79 // ---------------------------------------------------------------------
80 // =====================================================================
81 
82 int iProjectorJoseph::ReadConfigurationFile(const string& a_configurationFile)
83 {
84  // No options for joseph
85  ;
86  // Normal end
87  return 0;
88 }
89 
90 // =====================================================================
91 // ---------------------------------------------------------------------
92 // ---------------------------------------------------------------------
93 // =====================================================================
94 
95 int iProjectorJoseph::ReadOptionsList(const string& a_optionsList)
96 {
97  // No options for joseph
98  ;
99  // Normal end
100  return 0;
101 }
102 
103 // =====================================================================
104 // ---------------------------------------------------------------------
105 // ---------------------------------------------------------------------
106 // =====================================================================
107 
108 void iProjectorJoseph::ShowHelpSpecific()
109 {
110  cout << "This projector is a line projector that uses linear interpolation between pixels." << endl;
111  cout << "It is implemented from the following published paper:" << endl;
112  cout << "P. M. Joseph, \"An improved algorithm for reprojecting rays through pixel images\", IEEE Trans. Med. Imaging, vol. 1, pp. 192-6, 1982." << endl;
113 }
114 
115 // =====================================================================
116 // ---------------------------------------------------------------------
117 // ---------------------------------------------------------------------
118 // =====================================================================
119 
120 int iProjectorJoseph::CheckSpecificParameters()
121 {
122  // Nothing to check for this projector
123  ;
124  // Normal end
125  return 0;
126 }
127 
128 // =====================================================================
129 // ---------------------------------------------------------------------
130 // ---------------------------------------------------------------------
131 // =====================================================================
132 
133 int iProjectorJoseph::InitializeSpecific()
134 {
135  // Verbose
136  if (m_verbose>=2) Cout("iProjectorJoseph::InitializeSpecific() -> Use Joseph projector" << endl);
137 
138  // Allocate and compute boundaries (grid through voxel centers)
139  mp_boundariesX = new double[ mp_nbVox[ 0 ] + 2 ];
140  for( INTNB i = 0; i < mp_nbVox[ 0 ] + 2; ++i ) mp_boundariesX[ i ] = -((double)(mp_halfFOV[ 0 ])) + ((double)(mp_sizeVox[ 0 ])) * ( ((double)i) - 0.5 );
141  mp_boundariesY = new double[ mp_nbVox[ 1 ] + 2 ];
142  for( INTNB i = 0; i < mp_nbVox[ 1 ] + 2; ++i ) mp_boundariesY[ i ] = -((double)(mp_halfFOV[ 1 ])) + ((double)(mp_sizeVox[ 1 ])) * ( ((double)i) - 0.5 );
143  mp_boundariesZ = new double[ mp_nbVox[ 2 ] + 2 ];
144  for( INTNB i = 0; i < mp_nbVox[ 2 ] + 2; ++i ) mp_boundariesZ[ i ] = -((double)(mp_halfFOV[ 2 ])) + ((double)(mp_sizeVox[ 2 ])) * ( ((double)i) - 0.5 );
145 
146  // Allocating the mask buffer for the padded image space
147  INTNB nElts = ( mp_nbVox[ 0 ] + 2 ) * ( mp_nbVox[ 1 ] + 2 ) * ( mp_nbVox[ 2 ] + 2 );
148  mp_maskPad = new uint8_t[ nElts ];
149  ::memset( mp_maskPad, 0, sizeof( uint8_t ) * nElts );
150  for( INTNB k = 1; k < mp_nbVox[ 2 ] + 1; ++k )
151  {
152  for( INTNB j = 1; j < mp_nbVox[ 1 ] + 1; ++j )
153  {
154  for( INTNB i = 1; i < mp_nbVox[ 0 ] + 1; ++i )
155  {
156  mp_maskPad[ i + j * ( ( mp_nbVox[ 0 ] + 2 ) ) + k * ( mp_nbVox[ 0 ] + 2 ) * ( mp_nbVox[ 1 ] + 2 ) ] = 1;
157  }
158  }
159  }
160 
161  // Set the tolerance with respect to voxel sizes in each dimensions
162  double tolerance_factor = 1.e-4;
163  m_toleranceX = ((double)(mp_sizeVox[0])) * tolerance_factor;
164  m_toleranceY = ((double)(mp_sizeVox[1])) * tolerance_factor;
165  m_toleranceZ = ((double)(mp_sizeVox[2])) * tolerance_factor;
166 
167  // Normal end
168  return 0;
169 }
170 
171 // =====================================================================
172 // ---------------------------------------------------------------------
173 // ---------------------------------------------------------------------
174 // =====================================================================
175 
176 INTNB iProjectorJoseph::EstimateMaxNumberOfVoxelsPerLine()
177 {
178  // Find the maximum number of voxels along a given dimension
179  INTNB max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxX();
180  if (mp_ImageDimensionsAndQuantification->GetNbVoxY()>max_nb_voxels_in_dimension) max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxY();
181  if (mp_ImageDimensionsAndQuantification->GetNbVoxZ()>max_nb_voxels_in_dimension) max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxZ();
182  // We should have at most 4 voxels in a given plane, so multiply by 4
183  max_nb_voxels_in_dimension *= 4;
184  // Return the value
185  return max_nb_voxels_in_dimension;
186 }
187 
188 // =====================================================================
189 // ---------------------------------------------------------------------
190 // ---------------------------------------------------------------------
191 // =====================================================================
192 
193 int iProjectorJoseph::ProjectWithoutTOF(int a_direction, oProjectionLine* ap_ProjectionLine )
194 {
195  #ifdef CASTOR_DEBUG
196  if (!m_initialized)
197  {
198  Cerr("***** iProjectorJoseph::ProjectWithoutTOF() -> Called while not initialized !" << endl);
199  Exit(EXIT_DEBUG);
200  }
201  #endif
202 
203  #ifdef CASTOR_VERBOSE
204  if (m_verbose>=10)
205  {
206  string direction = "";
207  if (a_direction==FORWARD) direction = "forward";
208  else direction = "backward";
209  Cout("iProjectorJoseph::Project without TOF -> Project line '" << ap_ProjectionLine << "' in " << direction << " direction" << endl);
210  }
211  #endif
212 
213  // Get event positions
214  FLTNB* event1 = ap_ProjectionLine->GetPosition1();
215  FLTNB* event2 = ap_ProjectionLine->GetPosition2();
216 
217  // Distance between point event1 and event2
218  double const dX21( event2[ 0 ] - event1[ 0 ] );
219  double const dY21( event2[ 1 ] - event1[ 1 ] );
220  double const dZ21( event2[ 2 ] - event1[ 2 ] );
221 
222  // Compute absolute values of distances, they are used twice or more
223  double const fabs_dX21 = ::fabs(dX21);
224  double const fabs_dY21 = ::fabs(dY21);
225  double const fabs_dZ21 = ::fabs(dZ21);
226 
227  // Size of the voxels
228  double const szImg_X( mp_sizeVox[ 0 ] );
229  double const szImg_Y( mp_sizeVox[ 1 ] );
230  double const szImg_Z( mp_sizeVox[ 2 ] );
231 
232  // Get the number of the planes
233  INTNB const nPlane_X( mp_nbVox[ 0 ] );
234  INTNB const nPlane_Y( mp_nbVox[ 1 ] );
235  INTNB const nPlane_Z( mp_nbVox[ 2 ] );
236  INTNB const nPlane_XY( nPlane_X * nPlane_Y );
237 
238  // Get the number of the padded planes
239  INTNB const nPlane_X_Pad( nPlane_X + 2 );
240  INTNB const nPlane_Y_Pad( nPlane_Y + 2 );
241  INTNB const nPlane_XY_Pad( nPlane_X_Pad * nPlane_Y_Pad );
242 
243  // Coordinates of the first and last planes
244  double const xPlane_0 = -mp_halfFOV[ 0 ];
245  double const yPlane_0 = -mp_halfFOV[ 1 ];
246  double const zPlane_0 = -mp_halfFOV[ 2 ];
247  double const xPlane_N = mp_halfFOV[ 0 ];
248  double const yPlane_N = mp_halfFOV[ 1 ];
249  double const zPlane_N = mp_halfFOV[ 2 ];
250 
251  // Parameters for the index
252  INTNB index_x( 0 );
253  INTNB index_y( 0 );
254  INTNB index_z( 0 );
255 
256  // Parameter for the intersection point in each axis (at voxel centers)
257  // X
258  double intersection_x( event1[ 0 ] < xPlane_0 ?
259  xPlane_0 + szImg_X / 2.0 : xPlane_N - szImg_X / 2.0 );
260 
261  // Y
262  double intersection_y( event1[ 1 ] < yPlane_0 ?
263  yPlane_0 + szImg_Y / 2.0 : yPlane_N - szImg_Y / 2.0 );
264 
265  // Z
266  double intersection_z( event1[ 2 ] < zPlane_0 ?
267  zPlane_0 + szImg_Z / 2.0 : zPlane_N - szImg_Z / 2.0 );
268 
269  // Min. Max. coordinate for histo.
270  double const xPlane_0_min = xPlane_0 + szImg_X / 2.0;
271  double const yPlane_0_min = yPlane_0 + szImg_Y / 2.0;
272  double const zPlane_0_min = zPlane_0 + szImg_Z / 2.0;
273 
274  // Min. Max security. Do not take the last slice
275  double const xPlane_0_min_safety = xPlane_0 - ( szImg_X / 2.0 );
276  double const yPlane_0_min_safety = yPlane_0 - ( szImg_Y / 2.0 );
277  double const zPlane_0_min_safety = zPlane_0 - ( szImg_Z / 2.0 );
278  double const xPlane_N_max_safety = xPlane_N + ( szImg_X / 2.0 );
279  double const yPlane_N_max_safety = yPlane_N + ( szImg_Y / 2.0 );
280  double const zPlane_N_max_safety = zPlane_N + ( szImg_Z / 2.0 );
281 
282  // Incrementation of the plane
283  double const increment[] = {
284  event1[ 0 ] < xPlane_0 ? szImg_X : -szImg_X,
285  event1[ 1 ] < yPlane_0 ? szImg_Y : -szImg_Y,
286  event1[ 2 ] < zPlane_0 ? szImg_Z : -szImg_Z
287  };
288 
289  INTNB const increment_index[] = {
290  event1[ 0 ] < xPlane_0 ? 1 : -1,
291  event1[ 1 ] < yPlane_0 ? 1 : -1,
292  event1[ 2 ] < zPlane_0 ? 1 : -1
293  };
294 
295  INTNB const first_index[] = {
296  event1[ 0 ] < xPlane_0 ? 0 : nPlane_X - 1,
297  event1[ 1 ] < yPlane_0 ? 0 : nPlane_Y - 1,
298  event1[ 2 ] < zPlane_0 ? 0 : nPlane_Z - 1
299  };
300 
301  // Values of interpolation
302  double x1( 0.0 ); double x2( 0.0 );
303  double y1( 0.0 ); double y2( 0.0 );
304  double z1( 0.0 ); double z2( 0.0 );
305 
306  // Find the principal direction
307  // X direction
308  if( fabs_dX21 >= fabs_dY21 && fabs_dX21 >= fabs_dZ21 )
309  {
310  // Compute the weight to normalize the line
311  double const weight( ((double)(ap_ProjectionLine->GetLength())) * szImg_X / fabs_dX21 );
312 
313  // Take the increment X
314  double const incrementX = increment[ 0 ];
315 
316  // Compute first intersection in Y and Z
317  double const p1X_IntersectionX1 = ( intersection_x - event1[ 0 ] );
318  intersection_y = ( p1X_IntersectionX1 * dY21 / dX21 ) + event1[ 1 ];
319  intersection_z = ( p1X_IntersectionX1 * dZ21 / dX21 ) + event1[ 2 ];
320 
321  // Compute the increment (next intersection minus current intersection)
322  double const incrementY =
323  ( ( intersection_x + incrementX - event1[ 0 ] ) * dY21 / dX21 )
324  + event1[ 1 ] - intersection_y;
325  double const incrementZ =
326  ( ( intersection_x + incrementX - event1[ 0 ] ) * dZ21 / dX21 )
327  + event1[ 2 ] - intersection_z;
328 
329  // Take first index for x
330  index_x = first_index[ 0 ];
331 
332  // Loop over the Y-Z planes
333  for( INTNB i = 0; i < nPlane_X; ++i )
334  {
335  // Check if intersection outside of the image space
336  if( intersection_y < ( yPlane_0_min_safety + m_toleranceY )
337  || intersection_y >= ( yPlane_N_max_safety - m_toleranceY )
338  || intersection_z < ( zPlane_0_min_safety + m_toleranceZ )
339  || intersection_z >= ( zPlane_N_max_safety - m_toleranceZ ) )
340  {
341  index_x += increment_index[ 0 ];
342  intersection_y += incrementY;
343  intersection_z += incrementZ;
344  continue;
345  }
346 
347  // Find the index in the image
348  index_y = ::floor( ( intersection_y - yPlane_0_min ) / szImg_Y );
349  index_z = ::floor( ( intersection_z - zPlane_0_min ) / szImg_Z );
350 
351  // Compute distance between intersection point and Y and Z boundaries
352  y2 = ( intersection_y - mp_boundariesY[ index_y + 1 ] ) / szImg_Y;
353  y1 = 1.0 - y2;
354  z2 = ( intersection_z - mp_boundariesZ[ index_z + 1 ] ) / szImg_Z;
355  z1 = 1.0 - z2;
356 
357  // Index in the image space
358  INTNB index_final = index_x + index_y * nPlane_X + index_z * nPlane_XY;
359 
360  // Index in the padded image space
361  INTNB index_final_pad = ( index_x + 1 ) + ( index_y + 1 ) * nPlane_X_Pad + ( index_z + 1 ) * nPlane_XY_Pad;
362  INTNB maskIdx1 = mp_maskPad[ index_final_pad ];
363  INTNB maskIdx2 = mp_maskPad[ index_final_pad + nPlane_X_Pad ];
364  INTNB maskIdx3 = mp_maskPad[ index_final_pad + nPlane_XY_Pad ];
365  INTNB maskIdx4 = mp_maskPad[ index_final_pad + nPlane_X_Pad + nPlane_XY_Pad ];
366 
367  ap_ProjectionLine->AddVoxel(a_direction, ( index_final ) * maskIdx1, ( y1 * z1 ) * weight * maskIdx1);
368  ap_ProjectionLine->AddVoxel(a_direction, ( index_final + nPlane_X ) * maskIdx2, ( y2 * z1 ) * weight * maskIdx2);
369  ap_ProjectionLine->AddVoxel(a_direction, ( index_final + nPlane_XY ) * maskIdx3, ( y1 * z2 ) * weight * maskIdx3);
370  ap_ProjectionLine->AddVoxel(a_direction, ( index_final + nPlane_X + nPlane_XY ) * maskIdx4, ( y2 * z2 ) * weight * maskIdx4);
371 
372  // Increment
373  index_x += increment_index[ 0 ];
374  intersection_y += incrementY;
375  intersection_z += incrementZ;
376  }
377  } // Y direction
378  else if( fabs_dY21 > fabs_dX21 && fabs_dY21 >= fabs_dZ21 )
379  {
380  // Compute the weight to normalize the line
381  double const weight( ((double)(ap_ProjectionLine->GetLength())) * szImg_Y / fabs_dY21 );
382 
383  // Take the increment Y
384  double const incrementY = increment[ 1 ];
385 
386  // Compute first intersection in X and Z
387  double const p1Y_IntersectionY1 = ( intersection_y - event1[ 1 ] );
388  intersection_x = ( p1Y_IntersectionY1 * dX21 / dY21 ) + event1[ 0 ];
389  intersection_z = ( p1Y_IntersectionY1 * dZ21 / dY21 ) + event1[ 2 ];
390 
391  // Compute the increment (next intersection minus current intersection)
392  double const incrementX =
393  ( ( intersection_y + incrementY - event1[ 1 ] ) * dX21 / dY21 )
394  + event1[ 0 ] - intersection_x;
395  double const incrementZ =
396  ( ( intersection_y + incrementY - event1[ 1 ] ) * dZ21 / dY21 )
397  + event1[ 2 ] - intersection_z;
398 
399  // Take first index for y
400  index_y = first_index[ 1 ];
401 
402  // Loop over the X-Z planes
403  for( INTNB i = 0; i < nPlane_Y; ++i )
404  {
405  // Check if intersection outside of the image space
406  if( intersection_x < ( xPlane_0_min_safety + m_toleranceX )
407  || intersection_x >= ( xPlane_N_max_safety - m_toleranceX )
408  || intersection_z < ( zPlane_0_min_safety + m_toleranceZ )
409  || intersection_z >= ( zPlane_N_max_safety - m_toleranceZ ) )
410  {
411  index_y += increment_index[ 1 ];
412  intersection_x += incrementX;
413  intersection_z += incrementZ;
414  continue;
415  }
416 
417  // Find the index in the image
418  index_x = ::floor( ( intersection_x - xPlane_0_min ) / szImg_X );
419  index_z = ::floor( ( intersection_z - zPlane_0_min ) / szImg_Z );
420 
421  // Compute distance between intersection point and Y and Z boundaries
422  x2 = ( intersection_x - mp_boundariesX[ index_x + 1 ] ) / szImg_X;
423  x1 = 1.0 - x2;
424  z2 = ( intersection_z - mp_boundariesZ[ index_z + 1 ] ) / szImg_Z;
425  z1 = 1.0 - z2;
426 
427  // Index in the image space
428  uint32_t index_final = index_x + index_y * nPlane_X + index_z * nPlane_XY;
429 
430  // Index in the padded image space
431  INTNB index_final_pad = ( index_x + 1 ) + ( index_y + 1 ) * nPlane_X_Pad + ( index_z + 1 ) * nPlane_XY_Pad;
432  INTNB maskIdx1 = mp_maskPad[ index_final_pad ];
433  INTNB maskIdx2 = mp_maskPad[ index_final_pad + 1 ];
434  INTNB maskIdx3 = mp_maskPad[ index_final_pad + nPlane_XY_Pad ];
435  INTNB maskIdx4 = mp_maskPad[ index_final_pad + 1 + nPlane_XY_Pad ];
436 
437  ap_ProjectionLine->AddVoxel(a_direction, ( index_final ) * maskIdx1, ( x1 * z1 ) * weight * maskIdx1);
438  ap_ProjectionLine->AddVoxel(a_direction, ( index_final + 1 ) * maskIdx2, ( x2 * z1 ) * weight * maskIdx2);
439  ap_ProjectionLine->AddVoxel(a_direction, ( index_final + nPlane_XY ) * maskIdx3, ( x1 * z2 ) * weight * maskIdx3);
440  ap_ProjectionLine->AddVoxel(a_direction, ( index_final + 1 + nPlane_XY ) * maskIdx4, ( x2 * z2 ) * weight * maskIdx4);
441 
442  // Increment
443  index_y += increment_index[ 1 ];
444  intersection_x += incrementX;
445  intersection_z += incrementZ;
446  }
447  }
448  else // Z direction
449  {
450  // Compute the weight to normalize the line
451  double const weight( ((double)(ap_ProjectionLine->GetLength())) * szImg_Z / fabs_dZ21 );
452 
453  // Take the increment Z
454  double const incrementZ = increment[ 2 ];
455 
456  // Compute first intersection in X and Y
457  double const p1Z_IntersectionZ1 = ( intersection_z - event1[ 2 ] );
458  intersection_y = ( p1Z_IntersectionZ1 * dY21 / dZ21 ) + event1[ 1 ];
459  intersection_x = ( p1Z_IntersectionZ1 * dX21 / dZ21 ) + event1[ 0 ];
460 
461  // Compute the increment (next intersection minus current intersection)
462  double const incrementY =
463  ( ( intersection_z + incrementZ - event1[ 2 ] ) * dY21 / dZ21 )
464  + event1[ 1 ] - intersection_y;
465  double const incrementX =
466  ( ( intersection_z + incrementZ - event1[ 2 ] ) * dX21 / dZ21 )
467  + event1[ 0 ] - intersection_x;
468 
469  // Take first index for z
470  index_z = first_index[ 2 ];
471 
472  // Loop over the Y-Z planes
473  for( INTNB i = 0; i < nPlane_Z; ++i )
474  {
475  // Check if intersection outside of the image space
476  if( intersection_y < ( yPlane_0_min_safety + m_toleranceY )
477  || intersection_y >= ( yPlane_N_max_safety - m_toleranceY )
478  || intersection_x < ( xPlane_0_min_safety + m_toleranceX )
479  || intersection_x >= ( xPlane_N_max_safety - m_toleranceX ) )
480  {
481  index_z += increment_index[ 2 ];
482  intersection_y += incrementY;
483  intersection_x += incrementX;
484  continue;
485  }
486 
487  // Find the index in the image
488  index_y = ::floor( ( intersection_y - yPlane_0_min ) / szImg_Y );
489  index_x = ::floor( ( intersection_x - xPlane_0_min ) / szImg_X );
490 
491  // Compute distance between intersection point and Y and Z boundaries
492  y2 = ( intersection_y - mp_boundariesY[ index_y + 1 ] ) / szImg_Y;
493  y1 = 1.0 - y2;
494  x2 = ( intersection_x - mp_boundariesX[ index_x + 1 ] ) / szImg_X;
495  x1 = 1.0 - x2;
496 
497  // Index in the image space
498  uint32_t index_final = index_x + index_y * nPlane_X + index_z * nPlane_XY;
499 
500  // Index in the padded image space
501  INTNB index_final_pad = ( index_x + 1 ) + ( index_y + 1 ) * nPlane_X_Pad + ( index_z + 1 ) * nPlane_XY_Pad;
502  INTNB maskIdx1 = mp_maskPad[ index_final_pad ];
503  INTNB maskIdx2 = mp_maskPad[ index_final_pad + 1 ];
504  INTNB maskIdx3 = mp_maskPad[ index_final_pad + nPlane_X_Pad ];
505  INTNB maskIdx4 = mp_maskPad[ index_final_pad + 1 + nPlane_X_Pad ];
506 
507  ap_ProjectionLine->AddVoxel(a_direction, ( index_final ) * maskIdx1, ( x1 * y1 ) * weight * maskIdx1);
508  ap_ProjectionLine->AddVoxel(a_direction, ( index_final + 1 ) * maskIdx2, ( x2 * y1 ) * weight * maskIdx2);
509  ap_ProjectionLine->AddVoxel(a_direction, ( index_final + nPlane_X ) * maskIdx3, ( x1 * y2 ) * weight * maskIdx3);
510  ap_ProjectionLine->AddVoxel(a_direction, ( index_final + nPlane_X + 1 ) * maskIdx4, ( x2 * y2 ) * weight * maskIdx4);
511 
512  // Increment
513  index_z += increment_index[ 2 ];
514  intersection_y += incrementY;
515  intersection_x += incrementX;
516  }
517  }
518  return 0;
519 }
520 
521 // =====================================================================
522 // ---------------------------------------------------------------------
523 // ---------------------------------------------------------------------
524 // =====================================================================
525 
526 
527 
528 
529 
530 int iProjectorJoseph::ProjectWithTOFPos(int a_direction, oProjectionLine* ap_ProjectionLine)
531 {
532  #ifdef CASTOR_DEBUG
533  if (!m_initialized)
534  {
535  Cerr("***** iProjectorJoseph::ProjectWithTOFPos() -> Called while not initialized !" << endl);
536  Exit(EXIT_DEBUG);
537  }
538  #endif
539 
540  #ifdef CASTOR_VERBOSE
541  if (m_verbose>=10)
542  {
543  string direction = "";
544  if (a_direction==FORWARD) direction = "forward";
545  else direction = "backward";
546  Cout("iProjectorJoseph::Project with TOF pos -> Project line '" << ap_ProjectionLine << "' in " << direction << " direction" << endl);
547  }
548  #endif
549 
550  // Get event positions
551  FLTNB* event1 = ap_ProjectionLine->GetPosition1();
552  FLTNB* event2 = ap_ProjectionLine->GetPosition2();
553 
554  // LOR length
555  double lor_length = ap_ProjectionLine->GetLength();
556 
557  // Get TOF info
558  FLTNB tof_resolution = ap_ProjectionLine->GetTOFResolution();
559  FLTNB tof_measurement = ap_ProjectionLine->GetTOFMeasurement();
560 
561  // TOF standard deviation and truncation
562  double tof_sigma = tof_resolution * SPEED_OF_LIGHT * 0.5 / TWO_SQRT_TWO_LN_2;
563  double tof_half_span = tof_sigma * m_TOFnbSigmas;
564 
565  // convert delta time into delta length
566  double tof_deltaL = tof_measurement * SPEED_OF_LIGHT * 0.5;
567 
568  // special case for aberrant TOF measurements which must belong to scattered or random counts
569  // return an empty line because we know that the count does not depict the emitting object that we want to reconstruct
570  if ( fabs(tof_deltaL) > lor_length * 0.5)
571  {
572  return 0;
573  }
574 
575  // distance between the first event and the center of the Gaussian distribution along the LOR
576  double lor_tof_center = lor_length * 0.5 + tof_deltaL;
577 
578  // coefficient for conversion from erf use to integral of actual TOF function (Gaussian with maximum=1)
579  double tof_norm_coef = tof_sigma * sqrt(2*M_PI);
580 
581  // Distance between point event1 and event2
582  double const dX21( event2[ 0 ] - event1[ 0 ] );
583  double const dY21( event2[ 1 ] - event1[ 1 ] );
584  double const dZ21( event2[ 2 ] - event1[ 2 ] );
585 
586  // Compute absolute values of distances, they are used twice or more
587  double const fabs_dX21 = ::fabs(dX21);
588  double const fabs_dY21 = ::fabs(dY21);
589  double const fabs_dZ21 = ::fabs(dZ21);
590 
591  // Size of the voxels
592  double const szImg_X( mp_sizeVox[ 0 ] );
593  double const szImg_Y( mp_sizeVox[ 1 ] );
594  double const szImg_Z( mp_sizeVox[ 2 ] );
595 
596  // Get the number of the planes
597  INTNB const nPlane_X( mp_nbVox[ 0 ] );
598  INTNB const nPlane_Y( mp_nbVox[ 1 ] );
599  INTNB const nPlane_Z( mp_nbVox[ 2 ] );
600  INTNB const nPlane_XY( nPlane_X * nPlane_Y );
601 
602  // Get the number of the padded planes
603  INTNB const nPlane_X_Pad( nPlane_X + 2 );
604  INTNB const nPlane_Y_Pad( nPlane_Y + 2 );
605  INTNB const nPlane_XY_Pad( nPlane_X_Pad * nPlane_Y_Pad );
606 
607  // Coordinates of the first and last planes (passing through voxel edges)
608  double const xPlane_0 = -mp_halfFOV[ 0 ];
609  double const yPlane_0 = -mp_halfFOV[ 1 ];
610  double const zPlane_0 = -mp_halfFOV[ 2 ];
611  double const xPlane_N = mp_halfFOV[ 0 ];
612  double const yPlane_N = mp_halfFOV[ 1 ];
613  double const zPlane_N = mp_halfFOV[ 2 ];
614 
615  // Parameters for the index
616  INTNB index_x( 0 );
617  INTNB index_y( 0 );
618  INTNB index_z( 0 );
619 
620  // coordinates of the center of the Gaussian distribution (signed distances manage the order of event1/event2 along each axis)
621  double tof_center[] = {0,0,0};
622  // coordinates of the lower edge of the first voxel falling inside the truncated Gaussian distribution
623  double tof_edge_low[] = {0,0,0};
624  // coordinate of the upper edge of the last voxel falling inside the truncated Gaussian distribution
625  double tof_edge_high[] = {0,0,0};
626  // index along each axis of the first voxel falling inside the truncated Gaussian distribution
627  INTNB tof_index_low[] = {0,0,0};
628  // index along each axis of the last voxel falling inside the truncated Gaussian distribution
629  INTNB tof_index_high[] = {0,0,0};
630 
631  // TOF X
632  tof_center[0] = event1[ 0 ] + lor_tof_center * dX21 / lor_length;
633 
634  tof_edge_low[0] = tof_center[0] - tof_half_span * fabs_dX21 / lor_length;
635  // make the coordinate match the lowest voxel edge
636  tof_index_low[0] = max( (INTNB)::floor( (tof_edge_low[0] - xPlane_0) / szImg_X ), 0);
637  // if TOF outside of field of view, return empty line
638  if (tof_index_low[0]>nPlane_X-1) return 0;
639  tof_edge_low[0] = (double)tof_index_low[0] * szImg_X + xPlane_0;
640 
641  tof_edge_high[0] = tof_center[0] + tof_half_span * fabs_dX21 / lor_length;
642  // make the coordinate match the highest voxel edge
643  tof_index_high[0] = min( (INTNB)::floor( (tof_edge_high[0] - xPlane_0) / szImg_X ), nPlane_X-1);
644  // if TOF outside of field of view, return empty line
645  if (tof_index_high[0]<0) return 0;
646  tof_edge_high[0] = (double)(tof_index_high[0]+1) * szImg_X + xPlane_0;
647 
648  // TOF Y
649  tof_center[1] = event1[ 1 ] + lor_tof_center * dY21 / lor_length;
650 
651  tof_edge_low[1] = tof_center[1] - tof_half_span * fabs_dY21 / lor_length;
652  tof_index_low[1] = max( (INTNB)::floor( (tof_edge_low[1] - yPlane_0) / szImg_Y ), 0);
653  if (tof_index_low[1]>nPlane_Y-1) return 0;
654  tof_edge_low[1] = (double)tof_index_low[1] * szImg_Y + yPlane_0;
655 
656  tof_edge_high[1] = tof_center[1] + tof_half_span * fabs_dY21 / lor_length;
657  tof_index_high[1] = min( (INTNB)::floor( (tof_edge_high[1] - yPlane_0) / szImg_Y ), nPlane_Y-1);
658  if (tof_index_high[1]<0) return 0;
659  tof_edge_high[1] = (double)(tof_index_high[1]+1) * szImg_Y + yPlane_0;
660 
661  // TOF Z
662  tof_center[2] = event1[ 2 ] + lor_tof_center * dZ21 / lor_length;
663 
664  tof_edge_low[2] = tof_center[2] - tof_half_span * fabs_dZ21 / lor_length;
665  tof_index_low[2] = max( (INTNB)::floor( (tof_edge_low[2] - zPlane_0) / szImg_Z ), 0);
666  if (tof_index_low[2]>nPlane_Z-1) return 0;
667  tof_edge_low[2] = (double)tof_index_low[2] * szImg_Z + zPlane_0;
668 
669  tof_edge_high[2] = tof_center[2] + tof_half_span * fabs_dZ21 / lor_length;
670  tof_index_high[2] = min( (INTNB)::floor( (tof_edge_high[2] - zPlane_0) / szImg_Z ), nPlane_Z-1);
671  if (tof_index_high[2]<0) return 0;
672  tof_edge_high[2] = (double)(tof_index_high[2]+1) * szImg_Z + zPlane_0;
673 
674  // Parameter for the intersection point in each axis (intersection with voxel centers)
675  // X
676  double intersection_x( event1[ 0 ] < xPlane_0 ?
677  tof_edge_low[0] + szImg_X / 2.0 : tof_edge_high[0] - szImg_X / 2.0 );
678 
679  // Y
680  double intersection_y( event1[ 1 ] < yPlane_0 ?
681  tof_edge_low[1] + szImg_Y / 2.0 : tof_edge_high[1] - szImg_Y / 2.0 );
682 
683  // Z
684  double intersection_z( event1[ 2 ] < zPlane_0 ?
685  tof_edge_low[2] + szImg_Z / 2.0 : tof_edge_high[2] - szImg_Z / 2.0 );
686 
687  // Min. Max. coordinate for histo.
688  double const xPlane_0_min = xPlane_0 + szImg_X / 2.0;
689  double const yPlane_0_min = yPlane_0 + szImg_Y / 2.0;
690  double const zPlane_0_min = zPlane_0 + szImg_Z / 2.0;
691 
692  // Min. Max security. Do not take the last slice
693  double const xPlane_0_min_safety = xPlane_0 - ( szImg_X / 2.0 );
694  double const yPlane_0_min_safety = yPlane_0 - ( szImg_Y / 2.0 );
695  double const zPlane_0_min_safety = zPlane_0 - ( szImg_Z / 2.0 );
696  double const xPlane_N_max_safety = xPlane_N + ( szImg_X / 2.0 );
697  double const yPlane_N_max_safety = yPlane_N + ( szImg_Y / 2.0 );
698  double const zPlane_N_max_safety = zPlane_N + ( szImg_Z / 2.0 );
699 
700  // Incrementation of the plane
701  double const increment[] = {
702  event1[ 0 ] < xPlane_0 ? szImg_X : -szImg_X,
703  event1[ 1 ] < yPlane_0 ? szImg_Y : -szImg_Y,
704  event1[ 2 ] < zPlane_0 ? szImg_Z : -szImg_Z
705  };
706 
707  INTNB const increment_index[] = {
708  event1[ 0 ] < xPlane_0 ? 1 : -1,
709  event1[ 1 ] < yPlane_0 ? 1 : -1,
710  event1[ 2 ] < zPlane_0 ? 1 : -1
711  };
712 
713  INTNB const first_index[] = {
714  event1[ 0 ] < xPlane_0 ? (INTNB)tof_index_low[0] : (INTNB)tof_index_high[0],
715  event1[ 1 ] < yPlane_0 ? (INTNB)tof_index_low[1] : (INTNB)tof_index_high[1],
716  event1[ 2 ] < zPlane_0 ? (INTNB)tof_index_low[2] : (INTNB)tof_index_high[2],
717  };
718 
719 
720  // Values of interpolation
721  double x1( 0.0 ); double x2( 0.0 );
722  double y1( 0.0 ); double y2( 0.0 );
723  double z1( 0.0 ); double z2( 0.0 );
724 
725  // Find the principal direction
726  // X direction
727  if( fabs_dX21 >= fabs_dY21 && fabs_dX21 >= fabs_dZ21 )
728  {
729  // Take the increment X
730  double const incrementX = increment[ 0 ];
731 
732  // Compute first intersection in Y and Z
733  double const p1X_IntersectionX1 = ( intersection_x - event1[ 0 ] );
734  intersection_y = ( p1X_IntersectionX1 * dY21 / dX21 ) + event1[ 1 ];
735  intersection_z = ( p1X_IntersectionX1 * dZ21 / dX21 ) + event1[ 2 ];
736 
737  // Compute the increment (next intersection minus current intersection)
738  double const incrementY =
739  ( ( intersection_x + incrementX - event1[ 0 ] ) * dY21 / dX21 )
740  + event1[ 1 ] - intersection_y;
741  double const incrementZ =
742  ( ( intersection_x + incrementX - event1[ 0 ] ) * dZ21 / dX21 )
743  + event1[ 2 ] - intersection_z;
744 
745  // Take first index for x
746  index_x = first_index[ 0 ];
747 
748  // tof : erf of previous voxel edge - Gaussian center
749  double previous_edge_erf = erf( ( ( intersection_x - increment_index[ 0 ] * 0.5 *szImg_X - tof_center[0]) * lor_length / fabs_dX21 ) / ( sqrt(2.)*tof_sigma ) );
750  double next_edge_erf, tof_weight;
751 
752  // Loop over the Y-Z planes
753  for( INTNB i = 0; i < (INTNB)(tof_index_high[0] - tof_index_low[0] + 1); ++i )
754  {
755  // Check if intersection outside of the image space
756  if( intersection_y < ( yPlane_0_min_safety + m_toleranceY )
757  || intersection_y >= ( yPlane_N_max_safety - m_toleranceY )
758  || intersection_z < ( zPlane_0_min_safety + m_toleranceZ )
759  || intersection_z >= ( zPlane_N_max_safety - m_toleranceZ ) )
760  {
761  index_x += increment_index[ 0 ];
762  intersection_x += incrementX;
763  intersection_y += incrementY;
764  intersection_z += incrementZ;
765  continue;
766  }
767 
768  // Find the index in the image
769  index_y = ::floor( ( intersection_y - yPlane_0_min ) / szImg_Y );
770  index_z = ::floor( ( intersection_z - zPlane_0_min ) / szImg_Z );
771 
772  // Compute distance between intersection point and Y and Z boundaries
773  y2 = ( intersection_y - mp_boundariesY[ index_y + 1 ] ) / szImg_Y;
774  y1 = 1.0 - y2;
775  z2 = ( intersection_z - mp_boundariesZ[ index_z + 1 ] ) / szImg_Z;
776  z1 = 1.0 - z2;
777 
778  // Index in the image space
779  INTNB index_final = index_x + index_y * nPlane_X + index_z * nPlane_XY;
780 
781  // Index in the padded image space
782  INTNB index_final_pad = ( index_x + 1 ) + ( index_y + 1 ) * nPlane_X_Pad + ( index_z + 1 ) * nPlane_XY_Pad;
783  INTNB maskIdx1 = mp_maskPad[ index_final_pad ];
784  INTNB maskIdx2 = mp_maskPad[ index_final_pad + nPlane_X_Pad ];
785  INTNB maskIdx3 = mp_maskPad[ index_final_pad + nPlane_XY_Pad ];
786  INTNB maskIdx4 = mp_maskPad[ index_final_pad + nPlane_X_Pad + nPlane_XY_Pad ];
787 
788  // tof : erf of next voxel edge - Gaussian center
789  next_edge_erf = erf( ( ( intersection_x + increment_index[ 0 ] * 0.5 * szImg_X - tof_center[0] ) * lor_length / fabs_dX21 ) / (sqrt(2.)*tof_sigma) );
790  // integration of the Gaussian is done on the LOR portion matching the whole current voxel along X axis
791  tof_weight = 0.5 * fabs(previous_edge_erf - next_edge_erf) * tof_norm_coef ;
792  // keeping record of the previous edge, so as to save 1 erf computation
793  previous_edge_erf = next_edge_erf;
794 
795  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, 0, ( index_final ) * maskIdx1, ( y1 * z1 ) * tof_weight * maskIdx1);
796  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, 0, ( index_final + nPlane_X ) * maskIdx2, ( y2 * z1 ) * tof_weight * maskIdx2);
797  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, 0, ( index_final + nPlane_XY ) * maskIdx3, ( y1 * z2 ) * tof_weight * maskIdx3);
798  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, 0, ( index_final + nPlane_X + nPlane_XY ) * maskIdx4, ( y2 * z2 ) * tof_weight * maskIdx4);
799 
800  // Increment
801  index_x += increment_index[ 0 ];
802  intersection_x += incrementX;
803  intersection_y += incrementY;
804  intersection_z += incrementZ;
805  }
806  } // Y direction
807  else if( fabs_dY21 > fabs_dX21 && fabs_dY21 >= fabs_dZ21 )
808  {
809  // Take the increment Y
810  double const incrementY = increment[ 1 ];
811 
812  // Compute first intersection in X and Z
813  double const p1Y_IntersectionY1 = ( intersection_y - event1[ 1 ] );
814  intersection_x = ( p1Y_IntersectionY1 * dX21 / dY21 ) + event1[ 0 ];
815  intersection_z = ( p1Y_IntersectionY1 * dZ21 / dY21 ) + event1[ 2 ];
816 
817  // Compute the increment (next intersection minus current intersection)
818  double const incrementX =
819  ( ( intersection_y + incrementY - event1[ 1 ] ) * dX21 / dY21 )
820  + event1[ 0 ] - intersection_x;
821  double const incrementZ =
822  ( ( intersection_y + incrementY - event1[ 1 ] ) * dZ21 / dY21 )
823  + event1[ 2 ] - intersection_z;
824 
825  // Take first index for y
826  index_y = first_index[ 1 ];
827 
828  // tof : erf of previous voxel edge - Gaussian center
829  double previous_edge_erf = erf( ( ( intersection_y - increment_index[ 1 ] * 0.5 *szImg_Y - tof_center[1] ) * lor_length / fabs_dY21) / ( sqrt(2.)*tof_sigma ) );
830  double next_edge_erf, tof_weight;
831 
832  // Loop over the X-Z planes
833  for( INTNB i = 0; i < (INTNB)(tof_index_high[1] - tof_index_low[1] + 1); ++i )
834  {
835  // Check if intersection outside of the image space
836  if( intersection_x < ( xPlane_0_min_safety + m_toleranceX )
837  || intersection_x >= ( xPlane_N_max_safety - m_toleranceX )
838  || intersection_z < ( zPlane_0_min_safety + m_toleranceZ )
839  || intersection_z >= ( zPlane_N_max_safety - m_toleranceZ ) )
840  {
841  index_y += increment_index[ 1 ];
842  intersection_y += incrementY;
843  intersection_x += incrementX;
844  intersection_z += incrementZ;
845  continue;
846  }
847 
848  // Find the index in the image
849  index_x = ::floor( ( intersection_x - xPlane_0_min ) / szImg_X );
850  index_z = ::floor( ( intersection_z - zPlane_0_min ) / szImg_Z );
851 
852  // Compute distance between intersection point and Y and Z boundaries
853  x2 = ( intersection_x - mp_boundariesX[ index_x + 1 ] ) / szImg_X;
854  x1 = 1.0 - x2;
855  z2 = ( intersection_z - mp_boundariesZ[ index_z + 1 ] ) / szImg_Z;
856  z1 = 1.0 - z2;
857 
858  // Index in the image space
859  uint32_t index_final = index_x + index_y * nPlane_X + index_z * nPlane_XY;
860 
861  // Index in the padded image space
862  INTNB index_final_pad = ( index_x + 1 ) + ( index_y + 1 ) * nPlane_X_Pad + ( index_z + 1 ) * nPlane_XY_Pad;
863  INTNB maskIdx1 = mp_maskPad[ index_final_pad ];
864  INTNB maskIdx2 = mp_maskPad[ index_final_pad + 1 ];
865  INTNB maskIdx3 = mp_maskPad[ index_final_pad + nPlane_XY_Pad ];
866  INTNB maskIdx4 = mp_maskPad[ index_final_pad + 1 + nPlane_XY_Pad ];
867 
868  // tof : erf of next voxel edge - Gaussian center
869  next_edge_erf = erf( ( ( intersection_y + increment_index[ 1 ] * 0.5 * szImg_Y - tof_center[1] ) * lor_length / fabs_dY21 ) / (sqrt(2.)*tof_sigma) );
870  // integration of the Gaussian is done on the LOR portion matching the whole current voxel along X axis
871  tof_weight = 0.5 * fabs(previous_edge_erf - next_edge_erf) * tof_norm_coef;
872  // keeping record of the previous edge, so as to save 1 erf computation
873  previous_edge_erf = next_edge_erf;
874 
875  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, 0, ( index_final ) * maskIdx1, ( x1 * z1 ) * tof_weight * maskIdx1);
876  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, 0, ( index_final + 1 ) * maskIdx2, ( x2 * z1 ) * tof_weight * maskIdx2);
877  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, 0, ( index_final + nPlane_XY ) * maskIdx3, ( x1 * z2 ) * tof_weight * maskIdx3);
878  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, 0, ( index_final + 1 + nPlane_XY ) * maskIdx4, ( x2 * z2 ) * tof_weight * maskIdx4);
879 
880  // Increment
881  index_y += increment_index[ 1 ];
882  intersection_y += incrementY;
883  intersection_x += incrementX;
884  intersection_z += incrementZ;
885  }
886  }
887  else // Z direction
888  {
889  // Take the increment Z
890  double const incrementZ = increment[ 2 ];
891 
892  // Compute first intersection in X and Y
893  double const p1Z_IntersectionZ1 = ( intersection_z - event1[ 2 ] );
894  intersection_y = ( p1Z_IntersectionZ1 * dY21 / dZ21 ) + event1[ 1 ];
895  intersection_x = ( p1Z_IntersectionZ1 * dX21 / dZ21 ) + event1[ 0 ];
896 
897  // Compute the increment (next intersection minus current intersection)
898  double const incrementY =
899  ( ( intersection_z + incrementZ - event1[ 2 ] ) * dY21 / dZ21 )
900  + event1[ 1 ] - intersection_y;
901  double const incrementX =
902  ( ( intersection_z + incrementZ - event1[ 2 ] ) * dX21 / dZ21 )
903  + event1[ 0 ] - intersection_x;
904 
905  // Take first index for z
906  index_z = first_index[ 2 ];
907 
908  // tof : erf of previous voxel edge - Gaussian center
909  double previous_edge_erf = erf( ( ( intersection_z - increment_index[ 2 ] * 0.5 *szImg_Z - tof_center[2] ) * lor_length / fabs_dZ21) / ( sqrt(2.)*tof_sigma ) );
910  double next_edge_erf, tof_weight;
911 
912  // Loop over the Y-Z planes
913  for( INTNB i = 0; i < (INTNB)(tof_index_high[2] - tof_index_low[2] + 1); ++i )
914  {
915  // Check if intersection outside of the image space
916  if( intersection_y < ( yPlane_0_min_safety + m_toleranceY )
917  || intersection_y >= ( yPlane_N_max_safety - m_toleranceY )
918  || intersection_x < ( xPlane_0_min_safety + m_toleranceX )
919  || intersection_x >= ( xPlane_N_max_safety - m_toleranceX ) )
920  {
921  index_z += increment_index[ 2 ];
922  intersection_z += incrementZ;
923  intersection_y += incrementY;
924  intersection_x += incrementX;
925  continue;
926  }
927 
928  // Find the index in the image
929  index_y = ::floor( ( intersection_y - yPlane_0_min ) / szImg_Y );
930  index_x = ::floor( ( intersection_x - xPlane_0_min ) / szImg_X );
931 
932  // Compute distance between intersection point and Y and Z boundaries
933  y2 = ( intersection_y - mp_boundariesY[ index_y + 1 ] ) / szImg_Y;
934  y1 = 1.0 - y2;
935  x2 = ( intersection_x - mp_boundariesX[ index_x + 1 ] ) / szImg_X;
936  x1 = 1.0 - x2;
937 
938  // Index in the image space
939  uint32_t index_final = index_x + index_y * nPlane_X + index_z * nPlane_XY;
940 
941  // Index in the padded image space
942  INTNB index_final_pad = ( index_x + 1 ) + ( index_y + 1 ) * nPlane_X_Pad + ( index_z + 1 ) * nPlane_XY_Pad;
943  INTNB maskIdx1 = mp_maskPad[ index_final_pad ];
944  INTNB maskIdx2 = mp_maskPad[ index_final_pad + 1 ];
945  INTNB maskIdx3 = mp_maskPad[ index_final_pad + nPlane_X_Pad ];
946  INTNB maskIdx4 = mp_maskPad[ index_final_pad + 1 + nPlane_X_Pad ];
947 
948  // tof : erf of next voxel edge - Gaussian center
949  next_edge_erf = erf( ( ( intersection_z + increment_index[ 2 ] * 0.5 * szImg_Z - tof_center[2] ) * lor_length / fabs_dZ21 ) / (sqrt(2.)*tof_sigma) );
950  // integration of the Gaussian is done on the LOR portion matching the whole current voxel along X axis
951  tof_weight = 0.5 * fabs(previous_edge_erf - next_edge_erf) * tof_norm_coef;
952  // keeping record of the previous edge, so as to save 1 erf computation
953  previous_edge_erf = next_edge_erf;
954 
955  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, 0, ( index_final ) * maskIdx1, ( x1 * y1 ) * tof_weight * maskIdx1);
956  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, 0, ( index_final + 1 ) * maskIdx2, ( x2 * y1 ) * tof_weight * maskIdx2);
957  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, 0, ( index_final + nPlane_X ) * maskIdx3, ( x1 * y2 ) * tof_weight * maskIdx3);
958  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, 0, ( index_final + nPlane_X + 1 ) * maskIdx4, ( x2 * y2 ) * tof_weight * maskIdx4);
959 
960  // Increment
961  index_z += increment_index[ 2 ];
962  intersection_z += incrementZ;
963  intersection_y += incrementY;
964  intersection_x += incrementX;
965  }
966  }
967 
968  return 0;
969 
970 
971 }
972 
973 // =====================================================================
974 // ---------------------------------------------------------------------
975 // ---------------------------------------------------------------------
976 // =====================================================================
977 
978 int iProjectorJoseph::ProjectWithTOFBin(int a_direction, oProjectionLine* ap_ProjectionLine)
979 {
980  Cerr("***** iProjectorJoseph::ProjectWithTOFBin() -> Not yet implemented !" << endl);
981  return 1;
982 /*
983  #ifdef CASTOR_DEBUG
984  if (!m_initialized)
985  {
986  Cerr("***** iProjectorJoseph::ProjectWithTOFPos() -> Called while not initialized !" << endl);
987  Exit(EXIT_DEBUG);
988  }
989  #endif
990 
991  #ifdef CASTOR_VERBOSE
992  if (m_verbose>=10)
993  {
994  string direction = "";
995  if (a_direction==FORWARD) direction = "forward";
996  else direction = "backward";
997  Cout("iProjectorJoseph::Project with TOF pos -> Project line '" << ap_ProjectionLine << "' in " << direction << " direction" << endl);
998  }
999  #endif
1000 
1001  // Get event positions
1002  FLTNB* event1 = ap_ProjectionLine->GetPosition1();
1003  FLTNB* event2 = ap_ProjectionLine->GetPosition2();
1004 
1005  // LOR length
1006  double lor_length = ap_ProjectionLine->GetLength();
1007 
1008  // Get TOF info
1009  double tof_resolution = ap_ProjectionLine->GetTOFResolution();
1010  double tof_bin_size = ap_ProjectionLine->GetTOFBinSize();
1011  INTNB tof_nb_bins = ap_ProjectionLine->GetNbTOFBins();
1012 
1013  double tof_sigma = tof_resolution * SPEED_OF_LIGHT * 0.5 / TWO_SQRT_TWO_LN_2;
1014  double tof_half_span = tof_sigma * m_TOFnbSigmas;
1015 
1016  // minimum and maximum TOF bins for this LOR
1017  INTNB tof_bin_index_high = (INTNB)::floor(lor_length * 0.5 / (tof_bin_size * SPEED_OF_LIGHT));
1018  INTNB tof_bin_index_low = - tof_bin_index_high;
1019  tof_bin_index_low += tof_nb_bins / 2;
1020  tof_bin_index_high += tof_nb_bins / 2;
1021 
1022  // convert delta time for the lowest bin into delta length
1023  double tof_deltaL = (tof_bin_index_low-tof_nb_bins/2) * tof_bin_size * SPEED_OF_LIGHT * 0.5;
1024 
1025  // distance between the first event1 and the center of the Gaussian distribution of the lowest TOF bin along the LOR
1026  double lor_tof_center = lor_length * 0.5 + tof_deltaL;
1027 
1028  // coefficient for conversion from erf use to integral of actual TOF function
1029  double tof_norm_coef = lor_length / (double)(tof_bin_index_high - tof_bin_index_low + 1);
1030 
1031  // Distance between point event1 and event2
1032  double const dX21( event2[ 0 ] - event1[ 0 ] );
1033  double const dY21( event2[ 1 ] - event1[ 1 ] );
1034  double const dZ21( event2[ 2 ] - event1[ 2 ] );
1035 
1036  // Compute absolute values of distances, they are used twice or more
1037  double const fabs_dX21 = ::fabs(dX21);
1038  double const fabs_dY21 = ::fabs(dY21);
1039  double const fabs_dZ21 = ::fabs(dZ21);
1040 
1041  // Size of the voxels
1042  double const szImg_X( mp_sizeVox[ 0 ] );
1043  double const szImg_Y( mp_sizeVox[ 1 ] );
1044  double const szImg_Z( mp_sizeVox[ 2 ] );
1045 
1046  // Get the number of the planes
1047  INTNB const nPlane_X( mp_nbVox[ 0 ] );
1048  INTNB const nPlane_Y( mp_nbVox[ 1 ] );
1049  INTNB const nPlane_Z( mp_nbVox[ 2 ] );
1050  INTNB const nPlane_XY( nPlane_X * nPlane_Y );
1051 
1052  // Get the number of the padded planes
1053  INTNB const nPlane_X_Pad( nPlane_X + 2 );
1054  INTNB const nPlane_Y_Pad( nPlane_Y + 2 );
1055  INTNB const nPlane_XY_Pad( nPlane_X_Pad * nPlane_Y_Pad );
1056 
1057  // Coordinates of the first and last planes (passing through voxel edges)
1058  double const xPlane_0 = -mp_halfFOV[ 0 ];
1059  double const yPlane_0 = -mp_halfFOV[ 1 ];
1060  double const zPlane_0 = -mp_halfFOV[ 2 ];
1061  double const xPlane_N = mp_halfFOV[ 0 ];
1062  double const yPlane_N = mp_halfFOV[ 1 ];
1063  double const zPlane_N = mp_halfFOV[ 2 ];
1064 
1065  // Parameters for the index
1066  INTNB index_x( 0 );
1067  INTNB index_y( 0 );
1068  INTNB index_z( 0 );
1069 
1070  // Min. Max. coordinate for histo.
1071  double const xPlane_0_min = xPlane_0 + szImg_X / 2.0;
1072  double const yPlane_0_min = yPlane_0 + szImg_Y / 2.0;
1073  double const zPlane_0_min = zPlane_0 + szImg_Z / 2.0;
1074 
1075  // Min. Max security. Do not take the last slice
1076  double const xPlane_0_min_safety = xPlane_0 - ( szImg_X / 2.0 );
1077  double const yPlane_0_min_safety = yPlane_0 - ( szImg_Y / 2.0 );
1078  double const zPlane_0_min_safety = zPlane_0 - ( szImg_Z / 2.0 );
1079  double const xPlane_N_max_safety = xPlane_N + ( szImg_X / 2.0 );
1080  double const yPlane_N_max_safety = yPlane_N + ( szImg_Y / 2.0 );
1081  double const zPlane_N_max_safety = zPlane_N + ( szImg_Z / 2.0 );
1082 
1083  // Incrementation of the plane
1084  double const increment[] = {
1085  event1[ 0 ] < xPlane_0 ? szImg_X : -szImg_X,
1086  event1[ 1 ] < yPlane_0 ? szImg_Y : -szImg_Y,
1087  event1[ 2 ] < zPlane_0 ? szImg_Z : -szImg_Z
1088  };
1089 
1090  INTNB const increment_index[] = {
1091  event1[ 0 ] < xPlane_0 ? 1 : -1,
1092  event1[ 1 ] < yPlane_0 ? 1 : -1,
1093  event1[ 2 ] < zPlane_0 ? 1 : -1
1094  };
1095 
1096  // Values of interpolation
1097  double x1( 0.0 ); double x2( 0.0 );
1098  double y1( 0.0 ); double y2( 0.0 );
1099  double z1( 0.0 ); double z2( 0.0 );
1100 
1101 for (INTNB tof_bin=tof_bin_index_low; tof_bin<tof_bin_index_high; tof_bin++)
1102 {
1103 
1104  // coordinates of the center of the Gaussian distribution (signed distances manage the order of event1/event2 along each axis)
1105  double tof_center[] = {0,0,0};
1106  // coordinates of the lower edge of the first voxel falling inside the truncated Gaussian distribution
1107  double tof_edge_low[] = {0,0,0};
1108  // coordinate of the upper edge of the last voxel falling inside the truncated Gaussian distribution
1109  double tof_edge_high[] = {0,0,0};
1110  // voxel index along each axis of the the first voxel falling inside the truncated Gaussian distribution
1111  INTNB tof_index_low[] = {0,0,0};
1112  // voxel index along each axis of the upper edge of the last voxel falling inside the truncated Gaussian distribution
1113  INTNB tof_index_high[] = {0,0,0};
1114 
1115  // TOF X
1116  tof_center[0] = event1[ 0 ] + lor_tof_center * dX21 / lor_length;
1117 
1118  tof_edge_low[0] = tof_center[0] - tof_half_span * fabs_dX21 / lor_length;
1119  // make the coordinate match the lowest voxel edge
1120  tof_index_low[0] = max( (INTNB)::floor( (tof_edge_low[0] - xPlane_0) / szImg_X ), 0);
1121  // if LOR outside of field of view, return empty line
1122  if (tof_index_low[0]>nPlane_X-1) return 0;
1123  tof_edge_low[0] = (double)tof_index_low[0] * szImg_X + xPlane_0;
1124 
1125  tof_edge_high[0] = tof_center[0] + tof_half_span * fabs_dX21 / lor_length;
1126  // make the coordinate match the highest voxel edge
1127  tof_index_high[0] = min( (INTNB)::floor( (tof_edge_high[0] - xPlane_0) / szImg_X ), nPlane_X-1);
1128  // if LOR outside of field of view, return empty line
1129  if (tof_index_high[0]<0) return 0;
1130  tof_edge_high[0] = (double)(tof_index_high[0]+1) * szImg_X + xPlane_0;
1131 
1132  // TOF Y
1133  tof_center[1] = event1[ 1 ] + lor_tof_center * dY21 / lor_length;
1134 
1135  tof_edge_low[1] = tof_center[1] - tof_half_span * fabs_dY21 / lor_length;
1136  tof_index_low[1] = max( (INTNB)::floor( (tof_edge_low[1] - yPlane_0) / szImg_Y ), 0);
1137  if (tof_index_low[1]>nPlane_Y-1) return 0;
1138  tof_edge_low[1] = (double)tof_index_low[1] * szImg_Y + yPlane_0;
1139 
1140  tof_edge_high[1] = tof_center[1] + tof_half_span * fabs_dY21 / lor_length;
1141  tof_index_high[1] = min( (INTNB)::floor( (tof_edge_high[1] - yPlane_0) / szImg_Y ), nPlane_Y-1);
1142  if (tof_index_high[1]<0) return 0;
1143  tof_edge_high[1] = (double)(tof_index_high[1]+1) * szImg_Y + yPlane_0;
1144 
1145  // TOF Z
1146  tof_center[2] = event1[ 2 ] + lor_tof_center * dZ21 / lor_length;
1147 
1148  tof_edge_low[2] = tof_center[2] - tof_half_span * fabs_dZ21 / lor_length;
1149  tof_index_low[2] = max( (INTNB)::floor( (tof_edge_low[2] - zPlane_0) / szImg_Z ), 0);
1150  if (tof_index_low[2]>nPlane_Z-1) return 0;
1151  tof_edge_low[2] = (double)tof_index_low[2] * szImg_Z + zPlane_0;
1152 
1153  tof_edge_high[2] = tof_center[2] + tof_half_span * fabs_dZ21 / lor_length;
1154  tof_index_high[2] = min( (INTNB)::floor( (tof_edge_high[2] - zPlane_0) / szImg_Z ), nPlane_Z-1);
1155  if (tof_index_high[2]<0) return 0;
1156  tof_edge_high[2] = (double)(tof_index_high[2]+1) * szImg_Z + zPlane_0;
1157 
1158  // Parameter for the intersection point in each axis (intersection with voxel centers)
1159  // X
1160  double intersection_x( event1[ 0 ] < xPlane_0 ?
1161  tof_edge_low[0] + szImg_X / 2.0 : tof_edge_high[0] - szImg_X / 2.0 );
1162 
1163  // Y
1164  double intersection_y( event1[ 1 ] < yPlane_0 ?
1165  tof_edge_low[1] + szImg_Y / 2.0 : tof_edge_high[1] - szImg_Y / 2.0 );
1166 
1167  // Z
1168  double intersection_z( event1[ 2 ] < zPlane_0 ?
1169  tof_edge_low[2] + szImg_Z / 2.0 : tof_edge_high[2] - szImg_Z / 2.0 );
1170 
1171 
1172  INTNB const first_index[] = {
1173  event1[ 0 ] < xPlane_0 ? (INTNB)tof_index_low[0] : (INTNB)tof_index_high[0],
1174  event1[ 1 ] < yPlane_0 ? (INTNB)tof_index_low[1] : (INTNB)tof_index_high[1],
1175  event1[ 2 ] < zPlane_0 ? (INTNB)tof_index_low[2] : (INTNB)tof_index_high[2],
1176  };
1177 
1178  // Find the principal direction
1179  // X direction
1180  if( fabs_dX21 >= fabs_dY21 && fabs_dX21 >= fabs_dZ21 )
1181  {
1182  // Take the increment X
1183  double const incrementX = increment[ 0 ];
1184 
1185  // Compute first intersection in Y and Z
1186  double const p1X_IntersectionX1 = ( intersection_x - event1[ 0 ] );
1187  intersection_y = ( p1X_IntersectionX1 * dY21 / dX21 ) + event1[ 1 ];
1188  intersection_z = ( p1X_IntersectionX1 * dZ21 / dX21 ) + event1[ 2 ];
1189 
1190  // Compute the increment (next intersection minus current intersection)
1191  double const incrementY =
1192  ( ( intersection_x + incrementX - event1[ 0 ] ) * dY21 / dX21 )
1193  + event1[ 1 ] - intersection_y;
1194  double const incrementZ =
1195  ( ( intersection_x + incrementX - event1[ 0 ] ) * dZ21 / dX21 )
1196  + event1[ 2 ] - intersection_z;
1197 
1198  // Take first index for x
1199  index_x = first_index[ 0 ];
1200 
1201  // tof : erf of previous voxel edge - Gaussian center
1202  double previous_edge_erf = erf( ( fabs( intersection_x - increment_index[ 0 ] * 0.5 *szImg_X - tof_center[0]) * lor_length / fabs_dX21 ) / ( sqrt(2.)*tof_sigma ) );
1203 
1204  // Loop over the Y-Z planes
1205  for( INTNB i = 0; i < (INTNB)(tof_index_high[0] - tof_index_low[0] + 1); ++i )
1206  {
1207  // Check if intersection outside of the image space
1208  if( intersection_y < ( yPlane_0_min_safety + m_toleranceY )
1209  || intersection_y >= ( yPlane_N_max_safety - m_toleranceY )
1210  || intersection_z < ( zPlane_0_min_safety + m_toleranceZ )
1211  || intersection_z >= ( zPlane_N_max_safety - m_toleranceZ ) )
1212  {
1213  index_x += increment_index[ 0 ];
1214  intersection_x += incrementX;
1215  intersection_y += incrementY;
1216  intersection_z += incrementZ;
1217  continue;
1218  }
1219 
1220  // Find the index in the image
1221  index_y = ::floor( ( intersection_y - yPlane_0_min ) / szImg_Y );
1222  index_z = ::floor( ( intersection_z - zPlane_0_min ) / szImg_Z );
1223 
1224  // Compute distance between intersection point and Y and Z boundaries
1225  y2 = ( intersection_y - mp_boundariesY[ index_y + 1 ] ) / szImg_Y;
1226  y1 = 1.0 - y2;
1227  z2 = ( intersection_z - mp_boundariesZ[ index_z + 1 ] ) / szImg_Z;
1228  z1 = 1.0 - z2;
1229 
1230  // Index in the image space
1231  INTNB index_final = index_x + index_y * nPlane_X + index_z * nPlane_XY;
1232 
1233  // Index in the padded image space
1234  INTNB index_final_pad = ( index_x + 1 ) + ( index_y + 1 ) * nPlane_X_Pad + ( index_z + 1 ) * nPlane_XY_Pad;
1235  INTNB maskIdx1 = mp_maskPad[ index_final_pad ];
1236  INTNB maskIdx2 = mp_maskPad[ index_final_pad + nPlane_X_Pad ];
1237  INTNB maskIdx3 = mp_maskPad[ index_final_pad + nPlane_XY_Pad ];
1238  INTNB maskIdx4 = mp_maskPad[ index_final_pad + nPlane_X_Pad + nPlane_XY_Pad ];
1239 
1240  // tof : erf of next voxel edge - Gaussian center
1241  double next_edge_erf = erf( ( fabs( intersection_x + increment_index[ 0 ] * 0.5 * szImg_X - tof_center[0] ) * lor_length / fabs_dX21 ) / (sqrt(2.)*tof_sigma) );
1242  // integration of the Gaussian is done on the LOR portion matching the whole current voxel along X axis
1243  double tof_weight = 0.5 * fabs(previous_edge_erf - next_edge_erf) * tof_norm_coef ;
1244  // keeping record of the previous edge, so as to save 1 erf computation
1245  previous_edge_erf = next_edge_erf;
1246 
1247  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, tof_bin, ( index_final ) * maskIdx1, ( y1 * z1 ) * tof_weight * maskIdx1);
1248  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, tof_bin, ( index_final + nPlane_X ) * maskIdx2, ( y2 * z1 ) * tof_weight * maskIdx2);
1249  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, tof_bin, ( index_final + nPlane_XY ) * maskIdx3, ( y1 * z2 ) * tof_weight * maskIdx3);
1250  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, tof_bin, ( index_final + nPlane_X + nPlane_XY ) * maskIdx4, ( y2 * z2 ) * tof_weight * maskIdx4);
1251 
1252  // Increment
1253  index_x += increment_index[ 0 ];
1254  intersection_x += incrementX;
1255  intersection_y += incrementY;
1256  intersection_z += incrementZ;
1257  }
1258  } // Y direction
1259  else if( fabs_dY21 > fabs_dX21 && fabs_dY21 >= fabs_dZ21 )
1260  {
1261  // Take the increment Y
1262  double const incrementY = increment[ 1 ];
1263 
1264  // Compute first intersection in X and Z
1265  double const p1Y_IntersectionY1 = ( intersection_y - event1[ 1 ] );
1266  intersection_x = ( p1Y_IntersectionY1 * dX21 / dY21 ) + event1[ 0 ];
1267  intersection_z = ( p1Y_IntersectionY1 * dZ21 / dY21 ) + event1[ 2 ];
1268 
1269  // Compute the increment (next intersection minus current intersection)
1270  double const incrementX =
1271  ( ( intersection_y + incrementY - event1[ 1 ] ) * dX21 / dY21 )
1272  + event1[ 0 ] - intersection_x;
1273  double const incrementZ =
1274  ( ( intersection_y + incrementY - event1[ 1 ] ) * dZ21 / dY21 )
1275  + event1[ 2 ] - intersection_z;
1276 
1277  // Take first index for y
1278  index_y = first_index[ 1 ];
1279 
1280  // tof : erf of previous voxel edge - Gaussian center
1281  double previous_edge_erf = erf( ( fabs( intersection_y - increment_index[ 1 ] * 0.5 *szImg_Y - tof_center[1] ) * lor_length / fabs_dY21) / ( sqrt(2.)*tof_sigma ) );
1282 
1283  // Loop over the X-Z planes
1284  for( INTNB i = 0; i < (INTNB)(tof_index_high[1] - tof_index_low[1] + 1); ++i )
1285  {
1286  // Check if intersection outside of the image space
1287  if( intersection_x < ( xPlane_0_min_safety + m_toleranceX )
1288  || intersection_x >= ( xPlane_N_max_safety - m_toleranceX )
1289  || intersection_z < ( zPlane_0_min_safety + m_toleranceZ )
1290  || intersection_z >= ( zPlane_N_max_safety - m_toleranceZ ) )
1291  {
1292  index_y += increment_index[ 1 ];
1293  intersection_y += incrementY;
1294  intersection_x += incrementX;
1295  intersection_z += incrementZ;
1296  continue;
1297  }
1298 
1299  // Find the index in the image
1300  index_x = ::floor( ( intersection_x - xPlane_0_min ) / szImg_X );
1301  index_z = ::floor( ( intersection_z - zPlane_0_min ) / szImg_Z );
1302 
1303  // Compute distance between intersection point and Y and Z boundaries
1304  x2 = ( intersection_x - mp_boundariesX[ index_x + 1 ] ) / szImg_X;
1305  x1 = 1.0 - x2;
1306  z2 = ( intersection_z - mp_boundariesZ[ index_z + 1 ] ) / szImg_Z;
1307  z1 = 1.0 - z2;
1308 
1309  // Index in the image space
1310  uint32_t index_final = index_x + index_y * nPlane_X + index_z * nPlane_XY;
1311 
1312  // Index in the padded image space
1313  INTNB index_final_pad = ( index_x + 1 ) + ( index_y + 1 ) * nPlane_X_Pad + ( index_z + 1 ) * nPlane_XY_Pad;
1314  INTNB maskIdx1 = mp_maskPad[ index_final_pad ];
1315  INTNB maskIdx2 = mp_maskPad[ index_final_pad + 1 ];
1316  INTNB maskIdx3 = mp_maskPad[ index_final_pad + nPlane_XY_Pad ];
1317  INTNB maskIdx4 = mp_maskPad[ index_final_pad + 1 + nPlane_XY_Pad ];
1318 
1319  // tof : erf of next voxel edge - Gaussian center
1320  double next_edge_erf = erf( ( fabs( intersection_y + increment_index[ 1 ] * 0.5 * szImg_Y - tof_center[1] ) * lor_length / fabs_dY21 ) / (sqrt(2.)*tof_sigma) );
1321  // integration of the Gaussian is done on the LOR portion matching the whole current voxel along X axis
1322  double tof_weight = 0.5 * fabs(previous_edge_erf - next_edge_erf) * tof_norm_coef;
1323  // keeping record of the previous edge, so as to save 1 erf computation
1324  previous_edge_erf = next_edge_erf;
1325 
1326  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, tof_bin, ( index_final ) * maskIdx1, ( x1 * z1 ) * tof_weight * maskIdx1);
1327  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, tof_bin, ( index_final + 1 ) * maskIdx2, ( x2 * z1 ) * tof_weight * maskIdx2);
1328  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, tof_bin, ( index_final + nPlane_XY ) * maskIdx3, ( x1 * z2 ) * tof_weight * maskIdx3);
1329  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, tof_bin, ( index_final + 1 + nPlane_XY ) * maskIdx4, ( x2 * z2 ) * tof_weight * maskIdx4);
1330 
1331  // Increment
1332  index_y += increment_index[ 1 ];
1333  intersection_y += incrementY;
1334  intersection_x += incrementX;
1335  intersection_z += incrementZ;
1336  }
1337  }
1338  else // Z direction
1339  {
1340 
1341  // Take the increment Z
1342  double const incrementZ = increment[ 2 ];
1343 
1344  // Compute first intersection in X and Y
1345  double const p1Z_IntersectionZ1 = ( intersection_z - event1[ 2 ] );
1346  intersection_y = ( p1Z_IntersectionZ1 * dY21 / dZ21 ) + event1[ 1 ];
1347  intersection_x = ( p1Z_IntersectionZ1 * dX21 / dZ21 ) + event1[ 0 ];
1348 
1349  // Compute the increment (next intersection minus current intersection)
1350  double const incrementY =
1351  ( ( intersection_z + incrementZ - event1[ 2 ] ) * dY21 / dZ21 )
1352  + event1[ 1 ] - intersection_y;
1353  double const incrementX =
1354  ( ( intersection_z + incrementZ - event1[ 2 ] ) * dX21 / dZ21 )
1355  + event1[ 0 ] - intersection_x;
1356 
1357  // Take first index for z
1358  index_z = first_index[ 2 ];
1359 
1360  // tof : erf of previous voxel edge - Gaussian center
1361  double previous_edge_erf = erf( ( fabs( intersection_z - increment_index[ 2 ] * 0.5 *szImg_Z - tof_center[2] ) * lor_length / fabs_dZ21) / ( sqrt(2.)*tof_sigma ) );
1362 
1363  // Loop over the Y-Z planes
1364  for( INTNB i = 0; i < (INTNB)(tof_index_high[2] - tof_index_low[2] + 1); ++i )
1365  {
1366  // Check if intersection outside of the image space
1367  if( intersection_y < ( yPlane_0_min_safety + m_toleranceY )
1368  || intersection_y >= ( yPlane_N_max_safety - m_toleranceY )
1369  || intersection_x < ( xPlane_0_min_safety + m_toleranceX )
1370  || intersection_x >= ( xPlane_N_max_safety - m_toleranceX ) )
1371  {
1372  index_z += increment_index[ 2 ];
1373  intersection_z += incrementZ;
1374  intersection_y += incrementY;
1375  intersection_x += incrementX;
1376  continue;
1377  }
1378 
1379  // Find the index in the image
1380  index_y = ::floor( ( intersection_y - yPlane_0_min ) / szImg_Y );
1381  index_x = ::floor( ( intersection_x - xPlane_0_min ) / szImg_X );
1382 
1383  // Compute distance between intersection point and Y and Z boundaries
1384  y2 = ( intersection_y - mp_boundariesY[ index_y + 1 ] ) / szImg_Y;
1385  y1 = 1.0 - y2;
1386  x2 = ( intersection_x - mp_boundariesX[ index_x + 1 ] ) / szImg_X;
1387  x1 = 1.0 - x2;
1388 
1389  // Index in the image space
1390  uint32_t index_final = index_x + index_y * nPlane_X + index_z * nPlane_XY;
1391 
1392  // Index in the padded image space
1393  INTNB index_final_pad = ( index_x + 1 ) + ( index_y + 1 ) * nPlane_X_Pad + ( index_z + 1 ) * nPlane_XY_Pad;
1394  INTNB maskIdx1 = mp_maskPad[ index_final_pad ];
1395  INTNB maskIdx2 = mp_maskPad[ index_final_pad + 1 ];
1396  INTNB maskIdx3 = mp_maskPad[ index_final_pad + nPlane_X_Pad ];
1397  INTNB maskIdx4 = mp_maskPad[ index_final_pad + 1 + nPlane_X_Pad ];
1398 
1399  // tof : erf of next voxel edge - Gaussian center
1400  double next_edge_erf = erf( ( fabs( intersection_z + increment_index[ 2 ] * 0.5 * szImg_Z - tof_center[2] ) * lor_length / fabs_dZ21 ) / (sqrt(2.)*tof_sigma) );
1401  // integration of the Gaussian is done on the LOR portion matching the whole current voxel along X axis
1402  double tof_weight = 0.5 * fabs(previous_edge_erf - next_edge_erf) * tof_norm_coef;
1403  // keeping record of the previous edge, so as to save 1 erf computation
1404  previous_edge_erf = next_edge_erf;
1405 
1406  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, tof_bin, ( index_final ) * maskIdx1, ( x1 * y1 ) * tof_weight * maskIdx1);
1407  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, tof_bin, ( index_final + 1 ) * maskIdx2, ( x2 * y1 ) * tof_weight * maskIdx2);
1408  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, tof_bin, ( index_final + nPlane_X ) * maskIdx3, ( x1 * y2 ) * tof_weight * maskIdx3);
1409  ap_ProjectionLine->AddVoxelInTOFBin(a_direction, tof_bin, ( index_final + nPlane_X + 1 ) * maskIdx4, ( x2 * y2 ) * tof_weight * maskIdx4);
1410 
1411  // Increment
1412  index_z += increment_index[ 2 ];
1413  intersection_z += incrementZ;
1414  intersection_y += incrementY;
1415  intersection_x += incrementX;
1416  }
1417  }
1418 
1419  // update TOF center along the LOR for the next TOF bin
1420  lor_tof_center += tof_bin_size;
1421 
1422 }
1423 
1424  return 0;
1425 */
1426 
1427 }
1428 
1429 // =====================================================================
1430 // ---------------------------------------------------------------------
1431 // ---------------------------------------------------------------------
1432 // =====================================================================
vProjector::m_compatibleWithSPECTAttenuationCorrection
bool m_compatibleWithSPECTAttenuationCorrection
Definition: vProjector.hh:317
vProjector::mp_sizeVox
FLTNB mp_sizeVox[3]
Definition: vProjector.hh:301
vProjector::mp_nbVox
INTNB mp_nbVox[3]
Definition: vProjector.hh:302
FLTNB
#define FLTNB
Definition: gVariables.hh:55
iProjectorJoseph::InitializeSpecific
int InitializeSpecific()
This function is used to initialize specific stuff to the child projector.
Definition: iProjectorJoseph.cc:133
iProjectorJoseph::ShowHelpSpecific
void ShowHelpSpecific()
A function used to show help about the child module.
Definition: iProjectorJoseph.cc:108
iProjectorJoseph::ProjectWithTOFPos
int ProjectWithTOFPos(int a_direction, oProjectionLine *ap_ProjectionLine)
A function to project with TOF continuous information.
Definition: iProjectorJoseph.cc:530
TWO_SQRT_TWO_LN_2
#define TWO_SQRT_TWO_LN_2
Definition: gVariables.hh:72
vProjector::mp_ImageDimensionsAndQuantification
oImageDimensionsAndQuantification * mp_ImageDimensionsAndQuantification
Definition: vProjector.hh:307
vProjector
This class is designed to generically described any on-the-fly projector.
Definition: vProjector.hh:54
iProjectorJoseph::ProjectWithTOFBin
int ProjectWithTOFBin(int a_direction, oProjectionLine *ap_ProjectionLine)
A function to project with TOF binned information.
Definition: iProjectorJoseph.cc:978
iProjectorJoseph::ReadConfigurationFile
int ReadConfigurationFile(const string &a_configurationFile)
A function used to read options from a configuration file.
Definition: iProjectorJoseph.cc:82
oProjectionLine::GetTOFMeasurement
FLTNB GetTOFMeasurement()
This function is used to get the TOF measurement.
Definition: oProjectionLine.hh:336
oProjectionLine::GetLength
FLTNB GetLength()
This function is used to get the length of the line.
Definition: oProjectionLine.hh:245
iProjectorJoseph::ReadOptionsList
int ReadOptionsList(const string &a_optionsList)
A function used to read options from a list of options.
Definition: iProjectorJoseph.cc:95
iProjectorJoseph::CheckSpecificParameters
int CheckSpecificParameters()
A private function used to check the parameters settings specific to the child projector.
Definition: iProjectorJoseph.cc:120
Exit
void Exit(int code)
Definition: sOutputManager.cc:34
iProjectorJoseph::~iProjectorJoseph
~iProjectorJoseph()
The destructor of iProjectorJoseph.
Definition: iProjectorJoseph.cc:50
oProjectionLine::GetTOFResolution
FLTNB GetTOFResolution()
This function is used to get the TOF resolution.
Definition: oProjectionLine.hh:329
Cerr
#define Cerr(MESSAGE)
Definition: sOutputManager.hh:90
oProjectionLine::GetPosition1
FLTNB * GetPosition1()
This function is used to get the pointer to the mp_position1 (3-values tab).
Definition: oProjectionLine.hh:259
vProjector::m_initialized
bool m_initialized
Definition: vProjector.hh:323
oProjectionLine::AddVoxel
void AddVoxel(int a_direction, INTNB a_voxelIndice, FLTNB a_voxelWeight)
This function is used to add a voxel contribution to the line, assuming TOF bin 0 (i...
Definition: oProjectionLine.cc:610
iProjectorJoseph::EstimateMaxNumberOfVoxelsPerLine
INTNB EstimateMaxNumberOfVoxelsPerLine()
This function is used to compute and provide an estimate of the maximum number of voxels that could c...
Definition: iProjectorJoseph.cc:176
iProjectorJoseph.hh
Declaration of class iProjectorJoseph.
vProjector::m_TOFnbSigmas
FLTNB m_TOFnbSigmas
Definition: vProjector.hh:312
INTNB
#define INTNB
Definition: gVariables.hh:64
oProjectionLine
This class is designed to manage and store system matrix elements associated to a vEvent...
Definition: oProjectionLine.hh:70
iProjectorJoseph::ProjectWithoutTOF
int ProjectWithoutTOF(int a_direction, oProjectionLine *ap_ProjectionLine)
A function to project without TOF.
Definition: iProjectorJoseph.cc:193
sOutputManager.hh
Declaration of class sOutputManager.
vProjector::mp_halfFOV
FLTNB mp_halfFOV[3]
Definition: vProjector.hh:304
EXIT_DEBUG
#define EXIT_DEBUG
Definition: gVariables.hh:69
SPEED_OF_LIGHT
#define SPEED_OF_LIGHT
Definition: gVariables.hh:75
oProjectionLine::GetPosition2
FLTNB * GetPosition2()
This function is used to get the pointer to the mp_position2 (3-values tab).
Definition: oProjectionLine.hh:266
oImageDimensionsAndQuantification::GetNbVoxX
INTNB GetNbVoxX()
Get the number of voxels along the X axis.
Definition: oImageDimensionsAndQuantification.hh:576
oImageDimensionsAndQuantification::GetNbVoxZ
INTNB GetNbVoxZ()
Get the number of voxels along the Z axis.
Definition: oImageDimensionsAndQuantification.hh:590
FORWARD
#define FORWARD
Definition: oProjectionLine.hh:55
Cout
#define Cout(MESSAGE)
Definition: sOutputManager.hh:78
iProjectorJoseph::iProjectorJoseph
iProjectorJoseph()
The constructor of iProjectorJoseph.
Definition: iProjectorJoseph.cc:28
oProjectionLine::AddVoxelInTOFBin
void AddVoxelInTOFBin(int a_direction, int a_TOFBin, INTNB a_voxelIndice, FLTNB a_voxelWeight)
This function is used to add a voxel contribution to the line and provided TOF bin.
Definition: oProjectionLine.cc:554
oImageDimensionsAndQuantification::GetNbVoxY
INTNB GetNbVoxY()
Get the number of voxels along the Y axis.
Definition: oImageDimensionsAndQuantification.hh:583
Generated on Wed Jul 19 2017 17:02:32 for CASToR by   doxygen 1.8.6