CASToR  1.1
Tomographic Reconstruction (PET/SPECT)
 All Classes Files Functions Variables Typedefs Macros Groups Pages
oProjectionLine.cc
Go to the documentation of this file.
1 
2 /*
3  Implementation of class oProjectionLine
4 
5  - separators: done
6  - doxygen: done
7  - default initialization: done
8  - CASTOR_DEBUG:
9  - CASTOR_VERBOSE:
10  - TODO: Implement projection functions with attenuation for SPECT: ForwardProjectWithAttenuation() and BackwardProjectWithAttenuation()
11 */
12 
19 #include "oProjectionLine.hh"
20 #include "sOutputManager.hh"
21 #include "vProjector.hh"
22 
23 // =====================================================================
24 // ---------------------------------------------------------------------
25 // ---------------------------------------------------------------------
26 // =====================================================================
27 
29 {
30  // Default all values
31  m_verbose = 0;
32  m_checked = false;
33  m_initialized = false;
34  m_bedOffset = 0.;
35  m_nbTOFBins = -1;
36  m_currentTOFBin = 0;
37  m_TOFResolution = -1.;
38  m_TOFMeasurement = 0.;
39  mp_POI1 = NULL;
40  mp_POI2 = NULL;
41  mp_POIResolution = NULL;
42  m_UseMatchedProjectors = false;
43  mp_ForwardProjector = NULL;
44  mp_BackwardProjector = NULL;
45  m_length = 0.;
46  m2p_allocatedNbVoxels = NULL;
47  m2p_currentNbVoxels = NULL;
48  m3p_voxelIndices = NULL;
49  m3p_voxelWeights = NULL;
52  mp_position1 = NULL;
53  mp_position2 = NULL;
54  mp_orientation1 = NULL;
55  mp_orientation2 = NULL;
56  mp_bufferPosition1 = NULL;
57  mp_bufferPosition2 = NULL;
58  mp_bufferOrientation1 = NULL;
59  mp_bufferOrientation2 = NULL;
60  m_threadNumber = -1;
62  m_index1 = -1;
63  m_index2 = -1;
64 }
65 
66 // =====================================================================
67 // ---------------------------------------------------------------------
68 // ---------------------------------------------------------------------
69 // =====================================================================
70 
72 {
73  // Go through the destructor only if the object was initialized
74  if (m_initialized)
75  {
76  // Delete voxels lists
77  for (int b=0; b<m_nbTOFBins; b++)
78  {
79  if (m3p_voxelIndices[FORWARD][b])
80  {
82  else delete m3p_voxelIndices[FORWARD][b];
83  }
84  if (m3p_voxelWeights[FORWARD][b])
85  {
87  else delete m3p_voxelWeights[FORWARD][b];
88  }
90  {
91  if (m3p_voxelIndices[BACKWARD][b])
92  {
94  else delete m3p_voxelIndices[BACKWARD][b];
95  }
96  if (m3p_voxelWeights[BACKWARD][b])
97  {
99  else delete m3p_voxelWeights[BACKWARD][b];
100  }
101  }
102  }
103  // Delete rest
105  if (m3p_voxelWeights[FORWARD]) delete[] m3p_voxelWeights[FORWARD];
107  {
109  if (m3p_voxelWeights[BACKWARD]) delete[] m3p_voxelWeights[BACKWARD];
110  }
112  if (m2p_currentNbVoxels[FORWARD]) delete m2p_currentNbVoxels[FORWARD];
114  {
116  if (m2p_currentNbVoxels[BACKWARD]) delete m2p_currentNbVoxels[BACKWARD];
117  }
118  }
119 }
120 
121 // =====================================================================
122 // ---------------------------------------------------------------------
123 // ---------------------------------------------------------------------
124 // =====================================================================
125 
127 {
128  // Check mandatory parameters
129  if (m_nbTOFBins<=0)
130  {
131  Cerr("***** oProjectionLine::CheckParameters() -> Forbidden number of TOF bins (" << m_nbTOFBins << ") !" << endl);
132  Exit(EXIT_FAILURE);
133  }
135  {
136  Cerr("***** oProjectionLine::CheckParameters() -> Computation strategy incorrectly set !" << endl);
137  return 1;
138  }
139  if (mp_POIResolution==NULL)
140  {
141  Cerr("***** oProjectionLine::CheckParameters() -> POI resolution not set !" << endl);
142  return 1;
143  }
145  {
146  Cerr("***** oProjectionLine::CheckParameters() -> oImageDimensionsAndQuantification not set !" << endl);
147  return 1;
148  }
149  if (m_threadNumber<0)
150  {
151  Cerr("***** oProjectionLine::CheckParameters() -> The thread number associated to this line is not set !" << endl);
152  return 1;
153  }
154  // End
155  m_checked = true;
156  return 0;
157 }
158 
159 // =====================================================================
160 // ---------------------------------------------------------------------
161 // ---------------------------------------------------------------------
162 // =====================================================================
163 
165 {
166  // Forbid initialization without check
167  if (!m_checked)
168  {
169  Cerr("***** oProjectionLine::Initialize() -> Must call CheckParameters() before Initialize() !" << endl);
170  return 1;
171  }
172 
173  // Verbose
174  if (m_verbose>=4) Cout("oProjectionLine::Initialize() -> Initialize this projection line for thread " << m_threadNumber << endl);
175 
176  // Default LOR length
177  m_length = 0.;
178 
179  // Allocate pointers that depend on forward or backward projectors
180  m2p_allocatedNbVoxels = new INTNB*[2];
181  m2p_currentNbVoxels = new INTNB*[2];
182  m3p_voxelIndices = new INTNB**[2];
183  m3p_voxelWeights = new FLTNB**[2];
184 
185  // Allocate number of voxels and lists for forward projector
190 
191  // Allocate number of voxels and lists for backward projector
193  {
194  // If one single projector, then simply link pointers
199  }
200  else
201  {
202  // If separated projectors, then allocate
207  }
208 
209  // --------------------------------------------------------------------------------------
210  // Image computation strategy
212  {
213  // This image computation strategy is not compatible with system matrix projectors
214  if (mp_ForwardProjector==NULL || mp_BackwardProjector==NULL)
215  {
216  Cerr("***** oProjectionLine::Initialize() -> Image computation strategy is not compatible with the use of system matrix projectors !" << endl);
217  return 1;
218  }
219  // Verbose
220  if (m_verbose>=5) Cout(" --> Choose the image computation strategy" << endl);
221  // The number of voxels is held fixed to the image dimensions
223  // Set numbers and allocate
224  for (int b=0; b<m_nbTOFBins; b++)
225  {
226  // Just to remember the size of the image
227  m2p_allocatedNbVoxels[FORWARD][b] = nb_voxels;
228  m2p_currentNbVoxels[FORWARD][b] = nb_voxels;
229  // Voxel indices are useless with this computation strategy
230  m3p_voxelIndices[FORWARD][b] = NULL;
231  // Voxel weights
232  m3p_voxelWeights[FORWARD][b] = new FLTNB[nb_voxels];
233  // Do the same for backward projector if needed
235  {
236  // Just to remember the size of the image
237  m2p_allocatedNbVoxels[BACKWARD][b] = nb_voxels;
238  m2p_currentNbVoxels[BACKWARD][b] = nb_voxels;
239  // Voxel indices are useless with this computation strategy
240  m3p_voxelIndices[BACKWARD][b] = NULL;
241  // Voxel weights
242  m3p_voxelWeights[BACKWARD][b] = new FLTNB[nb_voxels];
243  }
244  }
245  }
246  // --------------------------------------------------------------------------------------
247  // Fixed list computation strategy
249  {
250  // Verbose
251  if (m_verbose>=5) Cout(" --> Choose the fixed list computation strategy" << endl);
252  // Set numbers and allocate
253  for (int b=0; b<m_nbTOFBins; b++)
254  {
255  // For the system matrix case
256  if (mp_ForwardProjector==NULL)
257  {
258  // Verbose
259  if (m_verbose>=5) Cout(" --> System matrix for forward projection" << endl);
260  // The number of allocated voxels is fixed with this strategy
263  // Voxel indices set to NULL for the moment because the pointer will be copied during the execution
264  m3p_voxelIndices[FORWARD][b] = NULL;
265  // Voxel weights also set to NULL
266  m3p_voxelWeights[FORWARD][b] = NULL;
267  }
268  // For the projector case
269  else
270  {
271  // The number of allocated voxels is fixed with this strategy
274  // Verbose
275  if (m_verbose>=5) Cout(" --> Projector for forward projection using " << m2p_allocatedNbVoxels[FORWARD][b] << " allocated voxels" << endl);
276  // Voxel indices
278  // Voxel weights
279  m3p_voxelWeights[FORWARD][b] = new FLTNB[m2p_allocatedNbVoxels[FORWARD][b]];
280  }
281  // Do the same for backward projector if needed
283  {
284  // For the system matrix case
285  if (mp_BackwardProjector==NULL)
286  {
287  // Verbose
288  if (m_verbose>=5) Cout(" --> System matrix for backward projection" << endl);
289  // The number of allocated voxels is fixed with this strategy
292  // Voxel indices set to NULL for the moment because the pointer will be copied during the execution
293  m3p_voxelIndices[BACKWARD][b] = NULL;
294  // Voxel weights also set to NULL
295  m3p_voxelWeights[BACKWARD][b] = NULL;
296  }
297  // For the projector case
298  else
299  {
300  // The number of allocated voxels is fixed with this strategy
303  // Verbose
304  if (m_verbose>=5) Cout(" --> Projector for backward projection using " << m2p_allocatedNbVoxels[BACKWARD][b] << " allocated voxels" << endl);
305  // Voxel indices
307  // Voxel weights
308  m3p_voxelWeights[BACKWARD][b] = new FLTNB[m2p_allocatedNbVoxels[BACKWARD][b]];
309  }
310  }
311  }
312  }
313  // --------------------------------------------------------------------------------------
314  // Adaptative list computation strategy
316  {
317  // This image computation strategy is not compatible with system matrix projectors
318  if (mp_ForwardProjector==NULL || mp_BackwardProjector==NULL)
319  {
320  Cerr("***** oProjectionLine::Initialize() -> Adaptative list computation strategy is not compatible with the use of system matrix projectors !" << endl);
321  return 1;
322  }
323  // The number of voxels is set to the diagonal as a first guess
325  // Verbose
326  if (m_verbose>=5) Cout(" --> Choose the adaptative list computation strategy, starting with " << nb_voxels << " allocated voxels" << endl);
327  // Set numbers and allocate
328  for (int b=0; b<m_nbTOFBins; b++)
329  {
330  // The number of allocated voxels is fixed with this strategy
331  m2p_allocatedNbVoxels[FORWARD][b] = nb_voxels;
333  // Voxel indices
334  m3p_voxelIndices[FORWARD][b] = (INTNB*)malloc(nb_voxels*sizeof(INTNB));
335  // Voxel weights
336  m3p_voxelWeights[FORWARD][b] = (FLTNB*)malloc(nb_voxels*sizeof(FLTNB));
337  // Do the same for backward projector if needed
339  {
340  // The number of allocated voxels is fixed with this strategy
341  m2p_allocatedNbVoxels[BACKWARD][b] = nb_voxels;
343  // Voxel indices
344  m3p_voxelIndices[BACKWARD][b] = (INTNB*)malloc(nb_voxels*sizeof(INTNB));
345  // Voxel weights
346  m3p_voxelWeights[BACKWARD][b] = (FLTNB*)malloc(nb_voxels*sizeof(FLTNB));
347  }
348  }
349  }
350 
351  // Allocate position and orientation of end points (as well as buffers)
352  mp_position1 = (FLTNB*)malloc(3*sizeof(FLTNB));
353  mp_position2 = (FLTNB*)malloc(3*sizeof(FLTNB));
354  mp_orientation1 = (FLTNB*)malloc(3*sizeof(FLTNB));
355  mp_orientation2 = (FLTNB*)malloc(3*sizeof(FLTNB));
356  mp_bufferPosition1 = (FLTNB*)malloc(3*sizeof(FLTNB));
357  mp_bufferPosition2 = (FLTNB*)malloc(3*sizeof(FLTNB));
358  mp_bufferOrientation1 = (FLTNB*)malloc(3*sizeof(FLTNB));
359  mp_bufferOrientation2 = (FLTNB*)malloc(3*sizeof(FLTNB));
360 
361  // End
362  m_initialized = true;
363  return 0;
364 }
365 
366 // =====================================================================
367 // ---------------------------------------------------------------------
368 // ---------------------------------------------------------------------
369 // =====================================================================
370 
372 {
373  #ifdef CASTOR_DEBUG
374  if (!m_initialized)
375  {
376  Cerr("***** oProjectionLine::ComputeLineLength() -> Called while not initialized !" << endl);
377  Exit(EXIT_DEBUG);
378  }
379  #endif
380 
381  #ifdef CASTOR_VERBOSE
382  if (m_verbose>=10) Cout("oProjectionLine::ComputeLineLength() -> Compute length of the line" << endl);
383  #endif
384 
388 }
389 
390 // =====================================================================
391 // ---------------------------------------------------------------------
392 // ---------------------------------------------------------------------
393 // =====================================================================
394 
396 {
397  #ifdef CASTOR_DEBUG
398  if (!m_initialized)
399  {
400  Cerr("***** oProjectionLine::NotEmptyLine() -> Called while not initialized !" << endl);
401  Exit(EXIT_DEBUG);
402  }
403  #endif
404 
405  #ifdef CASTOR_VERBOSE
406  if (m_verbose>=10) Cout("oProjectionLine::NotEmptyLine() -> Look if line is empty" << endl);
407  #endif
408 
409  // Just look in all TOF bins whether there are some contributing voxels
410  for (int b=0; b<m_nbTOFBins; b++) if (m2p_currentNbVoxels[FORWARD][b]!=0) return true;
411  return false;
412 }
413 
414 // =====================================================================
415 // ---------------------------------------------------------------------
416 // ---------------------------------------------------------------------
417 // =====================================================================
418 
420 {
421  #ifdef CASTOR_DEBUG
422  if (!m_initialized)
423  {
424  Cerr("***** oProjectionLine::Reset() -> Called while not initialized !" << endl);
425  Exit(EXIT_DEBUG);
426  }
427  #endif
428 
429  #ifdef CASTOR_VERBOSE
430  if (m_verbose>=10) Cout("oProjectionLine::Reset() -> Reset buffers of the line" << endl);
431  #endif
432 
433  // Set the length to zero
434  m_length = 0.;
435  // --------------------------------------------------------------------------------------
436  // Image computation strategy
438  {
439  // Loop on the TOF bins
440  for (int b=0; b<m_nbTOFBins; b++)
441  {
442  // Set all voxel weights to zero
443  for (int v=0; v<m2p_allocatedNbVoxels[FORWARD][b]; v++) m3p_voxelWeights[FORWARD][b][v] = 0.;
444  if (!m_UseMatchedProjectors) for (int v=0; v<m2p_allocatedNbVoxels[BACKWARD][b]; v++) m3p_voxelWeights[BACKWARD][b][v] = 0.;
445  }
446  }
447  // --------------------------------------------------------------------------------------
448  // List computation strategy
451  {
452  // Loop on the TOF bins
453  for (int b=0; b<m_nbTOFBins; b++)
454  {
455  // Only reset the number of current voxels, no need to reset matrices
458  }
459  }
460 }
461 
462 // =====================================================================
463 // ---------------------------------------------------------------------
464 // ---------------------------------------------------------------------
465 // =====================================================================
466 
468 {
469  #ifdef CASTOR_DEBUG
470  if (!m_initialized)
471  {
472  Cerr("***** oProjectionLine::ApplyOffset() -> Called while not initialized !" << endl);
473  Exit(EXIT_DEBUG);
474  }
475  #endif
476 
477  #ifdef CASTOR_VERBOSE
478  if (m_verbose>=10) Cout("oProjectionLine::ApplyOffset() -> Apply the global offset to the line end points" << endl);
479  #endif
480 
481  // Offset position 1
485  // Offset position 2
489 }
490 
491 // =====================================================================
492 // ---------------------------------------------------------------------
493 // ---------------------------------------------------------------------
494 // =====================================================================
495 
497 {
498  #ifdef CASTOR_DEBUG
499  if (!m_initialized)
500  {
501  Cerr("***** oProjectionLine::ApplyBedOffset() -> Called while not initialized !" << endl);
502  Exit(EXIT_DEBUG);
503  }
504  #endif
505 
506  #ifdef CASTOR_VERBOSE
507  if (m_verbose>=10) Cout("oProjectionLine::ApplyBedOffset() -> Apply the bed position offset to the line end points" << endl);
508  #endif
509 
510  // Offset position 1
512  // Offset position 2
514 }
515 
516 // =====================================================================
517 // ---------------------------------------------------------------------
518 // ---------------------------------------------------------------------
519 // =====================================================================
520 
521 INTNB oProjectionLine::GetVoxelIndex(int a_direction, int a_TOFBin, INTNB a_voxelInLine)
522 {
523  #ifdef CASTOR_DEBUG
524  if (!m_initialized)
525  {
526  Cerr("***** oProjectionLine::GetVoxelIndex() -> Called while not initialized !" << endl);
527  Exit(EXIT_DEBUG);
528  }
529  #endif
530 
531  #ifdef CASTOR_VERBOSE
532  if (m_verbose>=12)
533  {
534  string direction = "";
535  if (a_direction==FORWARD) direction = "forward";
536  else direction = "backward";
537  Cout("oProjectionLine::GetVoxelIndex() -> Get voxel index of voxel number " << a_voxelInLine << " in TOF bin " << a_TOFBin << " of " << direction << " projector" << endl);
538  }
539  #endif
540 
541  // If image computation strategy, then simply return the voxel index itself
543  return a_voxelInLine;
544  // Otherwise, look up for it in the voxel indices table
545  else
546  return m3p_voxelIndices[a_direction][a_TOFBin][a_voxelInLine];
547 }
548 
549 // =====================================================================
550 // ---------------------------------------------------------------------
551 // ---------------------------------------------------------------------
552 // =====================================================================
553 
554 void oProjectionLine::AddVoxelInTOFBin(int a_direction, int a_TOFBin, INTNB a_voxelIndex, FLTNB a_voxelWeight)
555 {
556  #ifdef CASTOR_DEBUG
557  if (!m_initialized)
558  {
559  Cerr("***** oProjectionLine::AddVoxelInTOFBin() -> Called while not initialized !" << endl);
560  Exit(EXIT_DEBUG);
561  }
562  #endif
563 
564  #ifdef CASTOR_VERBOSE
565  if (m_verbose>=12)
566  {
567  string direction = "";
568  if (a_direction==FORWARD) direction = "forward";
569  else direction = "backward";
570  Cout("oProjectionLine::AddVoxelInTOFBin() -> Add voxel index " << a_voxelIndex << " of weight " << a_voxelWeight <<
571  " into TOF bin " << a_TOFBin << " of " << direction << " projector" << endl);
572  }
573  #endif
574 
575  // Different way of doing for each computation strategy
577  {
578  // Simply add the contribution to this voxel index
579  m3p_voxelWeights[a_direction][a_TOFBin][a_voxelIndex] += a_voxelWeight;
580  }
582  {
583  // No buffer overflow checks for this computation strategy
584  m3p_voxelIndices[a_direction][a_TOFBin][m2p_currentNbVoxels[a_direction][a_TOFBin]] = a_voxelIndex;
585  m3p_voxelWeights[a_direction][a_TOFBin][m2p_currentNbVoxels[a_direction][a_TOFBin]] = a_voxelWeight;
586  m2p_currentNbVoxels[a_direction][a_TOFBin]++;
587  }
589  {
590  // Check if the number of contributing voxels is equal to the allocated size
591  if (m2p_currentNbVoxels[a_direction][a_TOFBin]==m2p_allocatedNbVoxels[a_direction][a_TOFBin])
592  {
593  // We realloc for one more voxel
594  m2p_allocatedNbVoxels[a_direction][a_TOFBin]++;
595  m3p_voxelIndices[a_direction][a_TOFBin] = (INTNB*)realloc(m3p_voxelIndices[a_direction][a_TOFBin],m2p_allocatedNbVoxels[a_direction][a_TOFBin]*sizeof(INTNB));
596  m3p_voxelWeights[a_direction][a_TOFBin] = (FLTNB*)realloc(m3p_voxelWeights[a_direction][a_TOFBin],m2p_allocatedNbVoxels[a_direction][a_TOFBin]*sizeof(FLTNB));
597  }
598  // Then register this voxel contribution
599  m3p_voxelIndices[a_direction][a_TOFBin][m2p_currentNbVoxels[a_direction][a_TOFBin]] = a_voxelIndex;
600  m3p_voxelWeights[a_direction][a_TOFBin][m2p_currentNbVoxels[a_direction][a_TOFBin]] = a_voxelWeight;
601  m2p_currentNbVoxels[a_direction][a_TOFBin]++;
602  }
603 }
604 
605 // =====================================================================
606 // ---------------------------------------------------------------------
607 // ---------------------------------------------------------------------
608 // =====================================================================
609 
610 void oProjectionLine::AddVoxel(int a_direction, INTNB a_voxelIndex, FLTNB a_voxelWeight)
611 {
612  #ifdef CASTOR_DEBUG
613  if (!m_initialized)
614  {
615  Cerr("***** oProjectionLine::AddVoxel() -> Called while not initialized !" << endl);
616  Exit(EXIT_DEBUG);
617  }
618  #endif
619 
620  #ifdef CASTOR_VERBOSE
621  if (m_verbose>=12)
622  {
623  string direction = "";
624  if (a_direction==FORWARD) direction = "forward";
625  else direction = "backward";
626  Cout("oProjectionLine::AddVoxel() -> Add voxel index " << a_voxelIndex << " of weight " << a_voxelWeight <<
627  " in " << direction << " projector" << endl);
628  }
629  #endif
630 
631  // No TOF here
632  int no_tof_bin = 0;
633 
634  // Different way of doing for each computation strategy
636  {
637  // Simply add the contribution to this voxel index
638  m3p_voxelWeights[a_direction][no_tof_bin][a_voxelIndex] += a_voxelWeight;
639  }
641  {
642  // No buffer overflow checks for this computation strategy
643  m3p_voxelIndices[a_direction][no_tof_bin][m2p_currentNbVoxels[a_direction][no_tof_bin]] = a_voxelIndex;
644  m3p_voxelWeights[a_direction][no_tof_bin][m2p_currentNbVoxels[a_direction][no_tof_bin]] = a_voxelWeight;
645  m2p_currentNbVoxels[a_direction][no_tof_bin]++;
646  }
648  {
649  // Check if the number of contributing voxels is equal to the allocated size
650  if (m2p_currentNbVoxels[a_direction][no_tof_bin]==m2p_allocatedNbVoxels[a_direction][no_tof_bin])
651  {
652  // We realloc for one more voxel
653  m2p_allocatedNbVoxels[a_direction][no_tof_bin]++;
654  m3p_voxelIndices[a_direction][no_tof_bin] = (INTNB*)realloc(m3p_voxelIndices[a_direction][no_tof_bin],m2p_allocatedNbVoxels[a_direction][no_tof_bin]*sizeof(INTNB));
655  m3p_voxelWeights[a_direction][no_tof_bin] = (FLTNB*)realloc(m3p_voxelWeights[a_direction][no_tof_bin],m2p_allocatedNbVoxels[a_direction][no_tof_bin]*sizeof(FLTNB));
656  }
657  // Then register this voxel contribution
658  m3p_voxelIndices[a_direction][no_tof_bin][m2p_currentNbVoxels[a_direction][no_tof_bin]] = a_voxelIndex;
659  m3p_voxelWeights[a_direction][no_tof_bin][m2p_currentNbVoxels[a_direction][no_tof_bin]] = a_voxelWeight;
660  m2p_currentNbVoxels[a_direction][no_tof_bin]++;
661  }
662 }
663 
664 // =====================================================================
665 // ---------------------------------------------------------------------
666 // ---------------------------------------------------------------------
667 // =====================================================================
668 
670 {
671  // The forward projection value
672  FLTNB value = 0.;
673  // If image is NULL, then assume 1 everywhere
674  if (ap_image==NULL)
675  {
676  // Loop over voxels inside the current TOF bin
677  for (int v=0; v<m2p_currentNbVoxels[FORWARD][m_currentTOFBin]; v++)
679  }
680  // Otherwise, project the image
681  else
682  {
684  {
685  // Loop over the whole image
686  for (int v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++)
687  value += m3p_voxelWeights[FORWARD][m_currentTOFBin][v] * ap_image[v];
688  }
689  else
690  {
691  // Loop over voxels inside the current TOF bin
692  for (int v=0; v<m2p_currentNbVoxels[FORWARD][m_currentTOFBin]; v++)
694  * ap_image[ m3p_voxelIndices[FORWARD][m_currentTOFBin][v] ];
695  }
696  }
697  // Apply multiplicative correction term
699  // Return the value
700  return value;
701 }
702 
703 // =====================================================================
704 // ---------------------------------------------------------------------
705 // ---------------------------------------------------------------------
706 // =====================================================================
707 
709 {
710  #ifdef CASTOR_DEBUG
711  // This function cannot be used by definition when using the image computation strategy
713  {
714  Cerr("***** oProjectionLine::ForwardProjectWithSPECTAttenuation() -> Cannot be used with an image computation strategy over the projection line !" << endl);
715  Exit(EXIT_FAILURE); // Have to exit
716  }
717  // This function cannot be used with an incompatible projector
719  {
720  Cerr("***** oProjectionLine::ForwardProjectWithSPECTAttenuation() -> Cannot be used with an incompatible projector !" << endl);
721  Exit(EXIT_FAILURE); // Have to exit
722  }
723  // Note that normally, the incompatibilities are checked before launching the reconstruction
724  #endif
725 
726  // The forward projection value
727  FLTNB value = 0.;
728 
729  // If image is NULL, then assume 1 everywhere in the emission object
730  if (ap_image==NULL)
731  {
732  // Loop over the element in the line
733  for (int v=0; v<m2p_currentNbVoxels[FORWARD][m_currentTOFBin]; v++)
734  {
735  FLTNB atn_sum = 0.0;
736  // Compute the attenuation of an element belonging to the line
737  for (int w=v; w<m2p_currentNbVoxels[FORWARD][m_currentTOFBin]; w++)
738  atn_sum += ap_attenuation[ m3p_voxelIndices[FORWARD][m_currentTOFBin][w] ] * m3p_voxelWeights[FORWARD][m_currentTOFBin][w];
739  // Take the inverse exponential to save a division after that (and convert from cm-1 to mm-1)
740  atn_sum = std::exp( -atn_sum *0.1 );
741 
742  value += m3p_voxelWeights[FORWARD][m_currentTOFBin][v] * atn_sum;
743  }
744  }
745  // Otherwise, project the image
746  else
747  {
748  // Loop over the element in the line
749  for (int v=0; v<m2p_currentNbVoxels[FORWARD][m_currentTOFBin]; v++)
750  {
751 
752  FLTNB atn_sum = 0.0;
753  // If attenuation image exist, compute the attenuation of an element belonging to the line
754  if(ap_attenuation != NULL)
755  for (int w=v; w<m2p_currentNbVoxels[FORWARD][m_currentTOFBin]; w++)
756  atn_sum += ap_attenuation[ m3p_voxelIndices[FORWARD][m_currentTOFBin][w] ] * m3p_voxelWeights[FORWARD][m_currentTOFBin][w];
757 
758  // Take the inverse exponential to save a division after that (and convert from cm-1 to mm-1)
759  atn_sum = std::exp( -atn_sum *0.1 );
760 
761  // Update projection
762  value += m3p_voxelWeights[FORWARD][m_currentTOFBin][v] * ap_image[ m3p_voxelIndices[FORWARD][m_currentTOFBin][v] ] * atn_sum;
763  }
764  }
765  // Apply multiplicative correction term
767  // Return the value
768  return value;
769 }
770 
771 // =====================================================================
772 // ---------------------------------------------------------------------
773 // ---------------------------------------------------------------------
774 // =====================================================================
775 
777 {
778  // Apply multiplicative correction term
779  a_value /= m_multiplicativeCorrection;
780  // Image-based strategy
782  {
783  // Loop over the whole image
784  for (int v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++)
785  ap_image[v] += a_value * m3p_voxelWeights[BACKWARD][m_currentTOFBin][v];
786  }
787  // List-based strategies
788  else
789  {
790  // Loop over voxels inside the current TOF bin
791  for (int v=0; v<m2p_currentNbVoxels[BACKWARD][m_currentTOFBin]; v++)
793  }
794 }
795 
796 // =====================================================================
797 // ---------------------------------------------------------------------
798 // ---------------------------------------------------------------------
799 // =====================================================================
800 
802 {
803  #ifdef CASTOR_DEBUG
804  // This function cannot be used by definition when using the image computation strategy
806  {
807  Cerr("***** oProjectionLine::BackwardProjectWithSPECTAttenuation() -> Cannot be used with an image computation strategy over the projection line !" << endl);
808  Exit(EXIT_FAILURE); // Have to exit
809  }
810  // This function cannot be used with an incompatible projector
812  {
813  Cerr("***** oProjectionLine::BackwardProjectWithSPECTAttenuation() -> Cannot be used with an incompatible projector !" << endl);
814  Exit(EXIT_FAILURE); // Have to exit
815  }
816  // Note that normally, the incompatibilities are checked before launching the reconstruction
817  #endif
818 
819  // Apply multiplicative correction term
820  a_value /= m_multiplicativeCorrection;
821 
822  // Loop over the element in the line
823  for (int v=0; v<m2p_currentNbVoxels[BACKWARD][m_currentTOFBin]; v++)
824  {
825  FLTNB atn_sum = 0.0;
826 
827  // If attenuation image exists, compute the attenuation of an element belonging to the line
828  if(ap_attenuation != NULL)
829  for (int w=v; w<m2p_currentNbVoxels[BACKWARD][m_currentTOFBin]; w++)
830  atn_sum += ap_attenuation[ m3p_voxelIndices[BACKWARD][m_currentTOFBin][w] ] * m3p_voxelWeights[BACKWARD][m_currentTOFBin][w];
831  // Take the inverse exponential to save a division after that (and convert from cm-1 to mm-1)
832  atn_sum = std::exp( -atn_sum *0.1 );
833  // Update image
834  ap_image[ m3p_voxelIndices[BACKWARD][m_currentTOFBin][v] ] += m3p_voxelWeights[BACKWARD][m_currentTOFBin][v] * a_value * atn_sum;
835  }
836 }
837 
838 // =====================================================================
839 // ---------------------------------------------------------------------
840 // ---------------------------------------------------------------------
841 // =====================================================================
842 
844 {
845  // Will be the cumulative integral
846  FLTNB integral = 0.;
847  // Simply loop over all voxels weights and sum them
848  for (int v=0; v<m2p_currentNbVoxels[a_direction][m_currentTOFBin]; v++)
849  integral += m3p_voxelWeights[a_direction][m_currentTOFBin][v];
850  // Return the result
851  return integral;
852 }
853 
~oProjectionLine()
The destructor of oProjectionLine.
#define IMAGE_COMPUTATION_STRATEGY
FLTNB ForwardProjectWithSPECTAttenuation(FLTNB *ap_attenuation, FLTNB *ap_image=NULL)
Forward projects the provided image for the current TOF bin with an inner loop on the attenuation (fo...
#define FLTNB
Definition: gVariables.hh:55
FLTNB * mp_bufferOrientation1
vProjector * mp_ForwardProjector
#define ADAPTATIVE_LIST_COMPUTATION_STRATEGY
bool GetCompatibilityWithSPECTAttenuationCorrection()
Definition: vProjector.hh:293
void Exit(int code)
oProjectionLine()
The constructor of oProjectionLine.
void ComputeLineLength()
Simply compute and update the m_length using the associated mp_position1 and mp_position2.
#define FIXED_LIST_COMPUTATION_STRATEGY
INTNB GetVoxelIndex(int a_direction, int a_TOFBin, INTNB a_voxelInLine)
This function is used to get the contributing voxel index of the provided direction, TOF bin and voxel rank.
INTNB *** m3p_voxelIndices
Declaration of class oProjectionLine.
void BackwardProject(FLTNB *ap_image, FLTNB a_value)
Simply backward projects the provided value inside the provided image, for the current TOF bin...
#define Cerr(MESSAGE)
void ApplyOffset()
Apply the offset of oImageDimensionsAndQuantification to the mp_position1 and mp_position2.
FLTNB *** m3p_voxelWeights
Declaration of class vProjector.
void BackwardProjectWithSPECTAttenuation(FLTNB *ap_attenuation, FLTNB *ap_image, FLTNB a_value)
Backward project the provided value inside the provided image with an inner loop on the attenuation (...
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...
void Reset()
Reset length and all the voxel indices and weights tabs.
bool NotEmptyLine()
This function is used to know if the line contains any voxel contribution.
FLTNB m_multiplicativeCorrection
FLTNB GetOffsetY()
Get the image offset along the Y axis, in mm.
#define INTNB
Definition: gVariables.hh:64
INTNB GetNbVoxDiagonal()
Get an estimation of the number of voxels along the image diagonal.
FLTNB GetOffsetZ()
Get the image offset along the Z axis, in mm.
Declaration of class sOutputManager.
FLTNB ForwardProject(FLTNB *ap_image=NULL)
Simply forward projects the provided image if not null, or else 1, for the current TOF bin...
INTNB GetNbVoxXYZ()
Get the total number of voxels.
FLTNB * mp_bufferPosition1
#define EXIT_DEBUG
Definition: gVariables.hh:69
FLTNB * mp_bufferOrientation2
FLTNB GetOffsetX()
Get the image offset along the X axis, in mm.
oImageDimensionsAndQuantification * mp_ImageDimensionsAndQuantification
#define FORWARD
#define Cout(MESSAGE)
INTNB ** m2p_allocatedNbVoxels
virtual INTNB EstimateMaxNumberOfVoxelsPerLine()
This function is used to compute and provide an estimate of the maximum number of voxels that could c...
Definition: vProjector.cc:266
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.
int CheckParameters()
A function used to check the parameters settings.
int Initialize()
A function used to initialize a bunch of stuff after parameters have been checked.
FLTNB * mp_bufferPosition2
void ApplyBedOffset()
Apply the bed offset of m_bedOffset to the mp_position1 and mp_position2.
INTNB ** m2p_currentNbVoxels
vProjector * mp_BackwardProjector
#define BACKWARD
FLTNB ComputeLineIntegral(int a_direction)
It simply computes the sum of all voxels contributions following the provided direction.