![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class oProjectionLine 00004 00005 - separators: done 00006 - doxygen: done 00007 - default initialization: done 00008 - CASTOR_DEBUG: 00009 - CASTOR_VERBOSE: 00010 - TODO: Implement projection functions with attenuation for SPECT: ForwardProjectWithAttenuation() and BackwardProjectWithAttenuation() 00011 */ 00012 00019 #include "oProjectionLine.hh" 00020 #include "sOutputManager.hh" 00021 #include "vProjector.hh" 00022 00023 // ===================================================================== 00024 // --------------------------------------------------------------------- 00025 // --------------------------------------------------------------------- 00026 // ===================================================================== 00027 00028 oProjectionLine::oProjectionLine() 00029 { 00030 // Default all values 00031 m_verbose = 0; 00032 m_checked = false; 00033 m_initialized = false; 00034 m_nbTOFBins = -1; 00035 m_currentTOFBin = 0; 00036 m_TOFResolution = -1.; 00037 m_TOFMeasurement = 0.; 00038 mp_POI1 = NULL; 00039 mp_POI2 = NULL; 00040 mp_POIResolution = NULL; 00041 m_UseMatchedProjectors = false; 00042 mp_ForwardProjector = NULL; 00043 mp_BackwardProjector = NULL; 00044 m_length = 0.; 00045 m2p_allocatedNbVoxels = NULL; 00046 m2p_currentNbVoxels = NULL; 00047 m3p_voxelIndices = NULL; 00048 m3p_voxelWeights = NULL; 00049 mp_ImageDimensionsAndQuantification = NULL; 00050 m_computationStrategy = -1; 00051 mp_position1 = NULL; 00052 mp_position2 = NULL; 00053 mp_orientation1 = NULL; 00054 mp_orientation2 = NULL; 00055 mp_bufferPosition1 = NULL; 00056 mp_bufferPosition2 = NULL; 00057 mp_bufferOrientation1 = NULL; 00058 mp_bufferOrientation2 = NULL; 00059 m_threadNumber = -1; 00060 m_multiplicativeCorrection = 1.; 00061 m_index1 = -1; 00062 m_index2 = -1; 00063 } 00064 00065 // ===================================================================== 00066 // --------------------------------------------------------------------- 00067 // --------------------------------------------------------------------- 00068 // ===================================================================== 00069 00070 oProjectionLine::~oProjectionLine() 00071 { 00072 // Go through the destructor only if the object was initialized 00073 if (m_initialized) 00074 { 00075 // Delete voxels lists 00076 for (int b=0; b<m_nbTOFBins; b++) 00077 { 00078 if (m3p_voxelIndices[FORWARD][b]) 00079 { 00080 if (m_computationStrategy==ADAPTATIVE_LIST_COMPUTATION_STRATEGY) free(m3p_voxelIndices[FORWARD][b]); 00081 else delete m3p_voxelIndices[FORWARD][b]; 00082 } 00083 if (m3p_voxelWeights[FORWARD][b]) 00084 { 00085 if (m_computationStrategy==ADAPTATIVE_LIST_COMPUTATION_STRATEGY) free(m3p_voxelWeights[FORWARD][b]); 00086 else delete m3p_voxelWeights[FORWARD][b]; 00087 } 00088 if (!m_UseMatchedProjectors) 00089 { 00090 if (m3p_voxelIndices[BACKWARD][b]) 00091 { 00092 if (m_computationStrategy==ADAPTATIVE_LIST_COMPUTATION_STRATEGY) free(m3p_voxelIndices[BACKWARD][b]); 00093 else delete m3p_voxelIndices[BACKWARD][b]; 00094 } 00095 if (m3p_voxelWeights[BACKWARD][b]) 00096 { 00097 if (m_computationStrategy==ADAPTATIVE_LIST_COMPUTATION_STRATEGY) free(m3p_voxelWeights[BACKWARD][b]); 00098 else delete m3p_voxelWeights[BACKWARD][b]; 00099 } 00100 } 00101 } 00102 // Delete rest 00103 if (m3p_voxelIndices[FORWARD]) delete[] m3p_voxelIndices[FORWARD]; 00104 if (m3p_voxelWeights[FORWARD]) delete[] m3p_voxelWeights[FORWARD]; 00105 if (!m_UseMatchedProjectors) 00106 { 00107 if (m3p_voxelIndices[BACKWARD]) delete[] m3p_voxelIndices[BACKWARD]; 00108 if (m3p_voxelWeights[BACKWARD]) delete[] m3p_voxelWeights[BACKWARD]; 00109 } 00110 if (m2p_allocatedNbVoxels[FORWARD]) delete m2p_allocatedNbVoxels[FORWARD]; 00111 if (m2p_currentNbVoxels[FORWARD]) delete m2p_currentNbVoxels[FORWARD]; 00112 if (!m_UseMatchedProjectors) 00113 { 00114 if (m2p_allocatedNbVoxels[BACKWARD]) delete m2p_allocatedNbVoxels[BACKWARD]; 00115 if (m2p_currentNbVoxels[BACKWARD]) delete m2p_currentNbVoxels[BACKWARD]; 00116 } 00117 } 00118 } 00119 00120 // ===================================================================== 00121 // --------------------------------------------------------------------- 00122 // --------------------------------------------------------------------- 00123 // ===================================================================== 00124 00125 int oProjectionLine::CheckParameters() 00126 { 00127 // Check mandatory parameters 00128 if (m_nbTOFBins<=0) 00129 { 00130 Cerr("***** oProjectionLine::CheckParameters() -> Forbidden number of TOF bins (" << m_nbTOFBins << ") !" << endl); 00131 Exit(EXIT_FAILURE); 00132 } 00133 if (m_computationStrategy!=IMAGE_COMPUTATION_STRATEGY && m_computationStrategy!=FIXED_LIST_COMPUTATION_STRATEGY && m_computationStrategy!=ADAPTATIVE_LIST_COMPUTATION_STRATEGY) 00134 { 00135 Cerr("***** oProjectionLine::CheckParameters() -> Computation strategy incorrectly set !" << endl); 00136 return 1; 00137 } 00138 if (mp_POIResolution==NULL) 00139 { 00140 Cerr("***** oProjectionLine::CheckParameters() -> POI resolution not set !" << endl); 00141 return 1; 00142 } 00143 if (mp_ImageDimensionsAndQuantification==NULL) 00144 { 00145 Cerr("***** oProjectionLine::CheckParameters() -> oImageDimensionsAndQuantification not set !" << endl); 00146 return 1; 00147 } 00148 if (m_threadNumber<0) 00149 { 00150 Cerr("***** oProjectionLine::CheckParameters() -> The thread number associated to this line is not set !" << endl); 00151 return 1; 00152 } 00153 // End 00154 m_checked = true; 00155 return 0; 00156 } 00157 00158 // ===================================================================== 00159 // --------------------------------------------------------------------- 00160 // --------------------------------------------------------------------- 00161 // ===================================================================== 00162 00163 int oProjectionLine::Initialize() 00164 { 00165 // Forbid initialization without check 00166 if (!m_checked) 00167 { 00168 Cerr("***** oProjectionLine::Initialize() -> Must call CheckParameters() before Initialize() !" << endl); 00169 return 1; 00170 } 00171 00172 // Verbose 00173 if (m_verbose>=4) Cout("oProjectionLine::Initialize() -> Initialize this projection line for thread " << m_threadNumber << endl); 00174 00175 // Default LOR length 00176 m_length = 0.; 00177 00178 // Allocate pointers that depend on forward or backward projectors 00179 m2p_allocatedNbVoxels = new INTNB*[2]; 00180 m2p_currentNbVoxels = new INTNB*[2]; 00181 m3p_voxelIndices = new INTNB**[2]; 00182 m3p_voxelWeights = new FLTNB**[2]; 00183 00184 // Allocate number of voxels and lists for forward projector 00185 m2p_allocatedNbVoxels[FORWARD] = new INTNB[m_nbTOFBins]; 00186 m2p_currentNbVoxels[FORWARD] = new INTNB[m_nbTOFBins]; 00187 m3p_voxelIndices[FORWARD] = new INTNB*[m_nbTOFBins]; 00188 m3p_voxelWeights[FORWARD] = new FLTNB*[m_nbTOFBins]; 00189 00190 // Allocate number of voxels and lists for backward projector 00191 if (m_UseMatchedProjectors) 00192 { 00193 // If one single projector, then simply link pointers 00194 m2p_allocatedNbVoxels[BACKWARD] = m2p_allocatedNbVoxels[FORWARD]; 00195 m2p_currentNbVoxels[BACKWARD] = m2p_currentNbVoxels[FORWARD]; 00196 m3p_voxelIndices[BACKWARD] = m3p_voxelIndices[FORWARD]; 00197 m3p_voxelWeights[BACKWARD] = m3p_voxelWeights[FORWARD]; 00198 } 00199 else 00200 { 00201 // If separated projectors, then allocate 00202 m2p_allocatedNbVoxels[BACKWARD] = new INTNB[m_nbTOFBins]; 00203 m2p_currentNbVoxels[BACKWARD] = new INTNB[m_nbTOFBins]; 00204 m3p_voxelIndices[BACKWARD] = new INTNB*[m_nbTOFBins]; 00205 m3p_voxelWeights[BACKWARD] = new FLTNB*[m_nbTOFBins]; 00206 } 00207 00208 // -------------------------------------------------------------------------------------- 00209 // Image computation strategy 00210 if (m_computationStrategy==IMAGE_COMPUTATION_STRATEGY) 00211 { 00212 // This image computation strategy is not compatible with system matrix projectors 00213 if (mp_ForwardProjector==NULL || mp_BackwardProjector==NULL) 00214 { 00215 Cerr("***** oProjectionLine::Initialize() -> Image computation strategy is not compatible with the use of system matrix projectors !" << endl); 00216 return 1; 00217 } 00218 // Verbose 00219 if (m_verbose>=5) Cout(" --> Choose the image computation strategy" << endl); 00220 // The number of voxels is held fixed to the image dimensions 00221 INTNB nb_voxels = mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); 00222 // Set numbers and allocate 00223 for (int b=0; b<m_nbTOFBins; b++) 00224 { 00225 // Just to remember the size of the image 00226 m2p_allocatedNbVoxels[FORWARD][b] = nb_voxels; 00227 m2p_currentNbVoxels[FORWARD][b] = nb_voxels; 00228 // Voxel indices are useless with this computation strategy 00229 m3p_voxelIndices[FORWARD][b] = NULL; 00230 // Voxel weights 00231 m3p_voxelWeights[FORWARD][b] = new FLTNB[nb_voxels]; 00232 // Do the same for backward projector if needed 00233 if (!m_UseMatchedProjectors) 00234 { 00235 // Just to remember the size of the image 00236 m2p_allocatedNbVoxels[BACKWARD][b] = nb_voxels; 00237 m2p_currentNbVoxels[BACKWARD][b] = nb_voxels; 00238 // Voxel indices are useless with this computation strategy 00239 m3p_voxelIndices[BACKWARD][b] = NULL; 00240 // Voxel weights 00241 m3p_voxelWeights[BACKWARD][b] = new FLTNB[nb_voxels]; 00242 } 00243 } 00244 } 00245 // -------------------------------------------------------------------------------------- 00246 // Fixed list computation strategy 00247 else if (m_computationStrategy==FIXED_LIST_COMPUTATION_STRATEGY) 00248 { 00249 // Verbose 00250 if (m_verbose>=5) Cout(" --> Choose the fixed list computation strategy" << endl); 00251 // Set numbers and allocate 00252 for (int b=0; b<m_nbTOFBins; b++) 00253 { 00254 // For the system matrix case 00255 if (mp_ForwardProjector==NULL) 00256 { 00257 // Verbose 00258 if (m_verbose>=5) Cout(" --> System matrix for forward projection" << endl); 00259 // The number of allocated voxels is fixed with this strategy 00260 m2p_allocatedNbVoxels[FORWARD][b] = 0; 00261 m2p_currentNbVoxels[FORWARD][b] = 0; 00262 // Voxel indices set to NULL for the moment because the pointer will be copied during the execution 00263 m3p_voxelIndices[FORWARD][b] = NULL; 00264 // Voxel weights also set to NULL 00265 m3p_voxelWeights[FORWARD][b] = NULL; 00266 } 00267 // For the projector case 00268 else 00269 { 00270 // The number of allocated voxels is fixed with this strategy 00271 m2p_allocatedNbVoxels[FORWARD][b] = mp_ForwardProjector->EstimateMaxNumberOfVoxelsPerLine(); 00272 m2p_currentNbVoxels[FORWARD][b] = 0; 00273 // Verbose 00274 if (m_verbose>=5) Cout(" --> Projector for forward projection using " << m2p_allocatedNbVoxels[FORWARD][b] << " allocated voxels" << endl); 00275 // Voxel indices 00276 m3p_voxelIndices[FORWARD][b] = new INTNB[m2p_allocatedNbVoxels[FORWARD][b]]; 00277 // Voxel weights 00278 m3p_voxelWeights[FORWARD][b] = new FLTNB[m2p_allocatedNbVoxels[FORWARD][b]]; 00279 } 00280 // Do the same for backward projector if needed 00281 if (!m_UseMatchedProjectors) 00282 { 00283 // For the system matrix case 00284 if (mp_BackwardProjector==NULL) 00285 { 00286 // Verbose 00287 if (m_verbose>=5) Cout(" --> System matrix for backward projection" << endl); 00288 // The number of allocated voxels is fixed with this strategy 00289 m2p_allocatedNbVoxels[BACKWARD][b] = 0; 00290 m2p_currentNbVoxels[BACKWARD][b] = 0; 00291 // Voxel indices set to NULL for the moment because the pointer will be copied during the execution 00292 m3p_voxelIndices[BACKWARD][b] = NULL; 00293 // Voxel weights also set to NULL 00294 m3p_voxelWeights[BACKWARD][b] = NULL; 00295 } 00296 // For the projector case 00297 else 00298 { 00299 // The number of allocated voxels is fixed with this strategy 00300 m2p_allocatedNbVoxels[BACKWARD][b] = mp_BackwardProjector->EstimateMaxNumberOfVoxelsPerLine(); 00301 m2p_currentNbVoxels[BACKWARD][b] = 0; 00302 // Verbose 00303 if (m_verbose>=5) Cout(" --> Projector for backward projection using " << m2p_allocatedNbVoxels[BACKWARD][b] << " allocated voxels" << endl); 00304 // Voxel indices 00305 m3p_voxelIndices[BACKWARD][b] = new INTNB[m2p_allocatedNbVoxels[BACKWARD][b]]; 00306 // Voxel weights 00307 m3p_voxelWeights[BACKWARD][b] = new FLTNB[m2p_allocatedNbVoxels[BACKWARD][b]]; 00308 } 00309 } 00310 } 00311 } 00312 // -------------------------------------------------------------------------------------- 00313 // Adaptative list computation strategy 00314 else if (m_computationStrategy==ADAPTATIVE_LIST_COMPUTATION_STRATEGY) 00315 { 00316 // This image computation strategy is not compatible with system matrix projectors 00317 if (mp_ForwardProjector==NULL || mp_BackwardProjector==NULL) 00318 { 00319 Cerr("***** oProjectionLine::Initialize() -> Adaptative list computation strategy is not compatible with the use of system matrix projectors !" << endl); 00320 return 1; 00321 } 00322 // The number of voxels is set to the diagonal as a first guess 00323 INTNB nb_voxels = mp_ImageDimensionsAndQuantification->GetNbVoxDiagonal(); 00324 // Verbose 00325 if (m_verbose>=5) Cout(" --> Choose the adaptative list computation strategy, starting with " << nb_voxels << " allocated voxels" << endl); 00326 // Set numbers and allocate 00327 for (int b=0; b<m_nbTOFBins; b++) 00328 { 00329 // The number of allocated voxels is fixed with this strategy 00330 m2p_allocatedNbVoxels[FORWARD][b] = nb_voxels; 00331 m2p_currentNbVoxels[FORWARD][b] = 0; 00332 // Voxel indices 00333 m3p_voxelIndices[FORWARD][b] = (INTNB*)malloc(nb_voxels*sizeof(INTNB)); 00334 // Voxel weights 00335 m3p_voxelWeights[FORWARD][b] = (FLTNB*)malloc(nb_voxels*sizeof(FLTNB)); 00336 // Do the same for backward projector if needed 00337 if (!m_UseMatchedProjectors) 00338 { 00339 // The number of allocated voxels is fixed with this strategy 00340 m2p_allocatedNbVoxels[BACKWARD][b] = nb_voxels; 00341 m2p_currentNbVoxels[BACKWARD][b] = 0; 00342 // Voxel indices 00343 m3p_voxelIndices[BACKWARD][b] = (INTNB*)malloc(nb_voxels*sizeof(INTNB)); 00344 // Voxel weights 00345 m3p_voxelWeights[BACKWARD][b] = (FLTNB*)malloc(nb_voxels*sizeof(FLTNB)); 00346 } 00347 } 00348 } 00349 00350 // Allocate position and orientation of end points (as well as buffers) 00351 mp_position1 = (FLTNB*)malloc(3*sizeof(FLTNB)); 00352 mp_position2 = (FLTNB*)malloc(3*sizeof(FLTNB)); 00353 mp_orientation1 = (FLTNB*)malloc(3*sizeof(FLTNB)); 00354 mp_orientation2 = (FLTNB*)malloc(3*sizeof(FLTNB)); 00355 mp_bufferPosition1 = (FLTNB*)malloc(3*sizeof(FLTNB)); 00356 mp_bufferPosition2 = (FLTNB*)malloc(3*sizeof(FLTNB)); 00357 mp_bufferOrientation1 = (FLTNB*)malloc(3*sizeof(FLTNB)); 00358 mp_bufferOrientation2 = (FLTNB*)malloc(3*sizeof(FLTNB)); 00359 00360 // End 00361 m_initialized = true; 00362 return 0; 00363 } 00364 00365 // ===================================================================== 00366 // --------------------------------------------------------------------- 00367 // --------------------------------------------------------------------- 00368 // ===================================================================== 00369 00370 void oProjectionLine::ComputeLineLength() 00371 { 00372 #ifdef CASTOR_DEBUG 00373 if (!m_initialized) 00374 { 00375 Cerr("***** oProjectionLine::ComputeLineLength() -> Called while not initialized !" << endl); 00376 Exit(EXIT_DEBUG); 00377 } 00378 #endif 00379 00380 #ifdef CASTOR_VERBOSE 00381 if (m_verbose>=10) Cout("oProjectionLine::ComputeLineLength() -> Compute length of the line" << endl); 00382 #endif 00383 00384 m_length = sqrt( (mp_position1[0]-mp_position2[0])*(mp_position1[0]-mp_position2[0]) + 00385 (mp_position1[1]-mp_position2[1])*(mp_position1[1]-mp_position2[1]) + 00386 (mp_position1[2]-mp_position2[2])*(mp_position1[2]-mp_position2[2]) ); 00387 } 00388 00389 // ===================================================================== 00390 // --------------------------------------------------------------------- 00391 // --------------------------------------------------------------------- 00392 // ===================================================================== 00393 00394 bool oProjectionLine::NotEmptyLine() 00395 { 00396 #ifdef CASTOR_DEBUG 00397 if (!m_initialized) 00398 { 00399 Cerr("***** oProjectionLine::NotEmptyLine() -> Called while not initialized !" << endl); 00400 Exit(EXIT_DEBUG); 00401 } 00402 #endif 00403 00404 #ifdef CASTOR_VERBOSE 00405 if (m_verbose>=10) Cout("oProjectionLine::NotEmptyLine() -> Look if line is empty" << endl); 00406 #endif 00407 00408 // Just look in all TOF bins whether there are some contributing voxels 00409 for (int b=0; b<m_nbTOFBins; b++) if (m2p_currentNbVoxels[FORWARD][b]!=0) return true; 00410 return false; 00411 } 00412 00413 // ===================================================================== 00414 // --------------------------------------------------------------------- 00415 // --------------------------------------------------------------------- 00416 // ===================================================================== 00417 00418 void oProjectionLine::Reset() 00419 { 00420 #ifdef CASTOR_DEBUG 00421 if (!m_initialized) 00422 { 00423 Cerr("***** oProjectionLine::Reset() -> Called while not initialized !" << endl); 00424 Exit(EXIT_DEBUG); 00425 } 00426 #endif 00427 00428 #ifdef CASTOR_VERBOSE 00429 if (m_verbose>=10) Cout("oProjectionLine::Reset() -> Reset buffers of the line" << endl); 00430 #endif 00431 00432 // Set the length to zero 00433 m_length = 0.; 00434 // -------------------------------------------------------------------------------------- 00435 // Image computation strategy 00436 if (m_computationStrategy==IMAGE_COMPUTATION_STRATEGY) 00437 { 00438 // Loop on the TOF bins 00439 for (int b=0; b<m_nbTOFBins; b++) 00440 { 00441 // Set all voxel weights to zero 00442 for (int v=0; v<m2p_allocatedNbVoxels[FORWARD][b]; v++) m3p_voxelWeights[FORWARD][b][v] = 0.; 00443 if (!m_UseMatchedProjectors) for (int v=0; v<m2p_allocatedNbVoxels[BACKWARD][b]; v++) m3p_voxelWeights[BACKWARD][b][v] = 0.; 00444 } 00445 } 00446 // -------------------------------------------------------------------------------------- 00447 // List computation strategy 00448 else if ( m_computationStrategy==FIXED_LIST_COMPUTATION_STRATEGY || 00449 m_computationStrategy==ADAPTATIVE_LIST_COMPUTATION_STRATEGY ) 00450 { 00451 // Loop on the TOF bins 00452 for (int b=0; b<m_nbTOFBins; b++) 00453 { 00454 // Only reset the number of current voxels, no need to reset matrices 00455 m2p_currentNbVoxels[FORWARD][b] = 0; 00456 if (!m_UseMatchedProjectors) m2p_currentNbVoxels[BACKWARD][b] = 0; 00457 } 00458 } 00459 } 00460 00461 // ===================================================================== 00462 // --------------------------------------------------------------------- 00463 // --------------------------------------------------------------------- 00464 // ===================================================================== 00465 00466 void oProjectionLine::ApplyOffset() 00467 { 00468 #ifdef CASTOR_DEBUG 00469 if (!m_initialized) 00470 { 00471 Cerr("***** oProjectionLine::ApplyOffset() -> Called while not initialized !" << endl); 00472 Exit(EXIT_DEBUG); 00473 } 00474 #endif 00475 00476 #ifdef CASTOR_VERBOSE 00477 if (m_verbose>=10) Cout("oProjectionLine::ApplyOffset() -> Apply the global offset to the line end points" << endl); 00478 #endif 00479 00480 // Offset position 1 00481 mp_position1[0] += mp_ImageDimensionsAndQuantification->GetOffsetX(); 00482 mp_position1[1] += mp_ImageDimensionsAndQuantification->GetOffsetY(); 00483 mp_position1[2] += mp_ImageDimensionsAndQuantification->GetOffsetZ(); 00484 // Offset position 2 00485 mp_position2[0] += mp_ImageDimensionsAndQuantification->GetOffsetX(); 00486 mp_position2[1] += mp_ImageDimensionsAndQuantification->GetOffsetY(); 00487 mp_position2[2] += mp_ImageDimensionsAndQuantification->GetOffsetZ(); 00488 } 00489 00490 // ===================================================================== 00491 // --------------------------------------------------------------------- 00492 // --------------------------------------------------------------------- 00493 // ===================================================================== 00494 00495 INTNB oProjectionLine::GetVoxelIndex(int a_direction, int a_TOFBin, INTNB a_voxelInLine) 00496 { 00497 #ifdef CASTOR_DEBUG 00498 if (!m_initialized) 00499 { 00500 Cerr("***** oProjectionLine::GetVoxelIndex() -> Called while not initialized !" << endl); 00501 Exit(EXIT_DEBUG); 00502 } 00503 #endif 00504 00505 #ifdef CASTOR_VERBOSE 00506 if (m_verbose>=12) 00507 { 00508 string direction = ""; 00509 if (a_direction==FORWARD) direction = "forward"; 00510 else direction = "backward"; 00511 Cout("oProjectionLine::GetVoxelIndex() -> Get voxel index of voxel number " << a_voxelInLine << " in TOF bin " << a_TOFBin << " of " << direction << " projector" << endl); 00512 } 00513 #endif 00514 00515 // If image computation strategy, then simply return the voxel index itself 00516 if (m_computationStrategy==IMAGE_COMPUTATION_STRATEGY) 00517 return a_voxelInLine; 00518 // Otherwise, look up for it in the voxel indices table 00519 else 00520 return m3p_voxelIndices[a_direction][a_TOFBin][a_voxelInLine]; 00521 } 00522 00523 // ===================================================================== 00524 // --------------------------------------------------------------------- 00525 // --------------------------------------------------------------------- 00526 // ===================================================================== 00527 00528 void oProjectionLine::AddVoxelInTOFBin(int a_direction, int a_TOFBin, INTNB a_voxelIndex, FLTNB a_voxelWeight) 00529 { 00530 #ifdef CASTOR_DEBUG 00531 if (!m_initialized) 00532 { 00533 Cerr("***** oProjectionLine::AddVoxelInTOFBin() -> Called while not initialized !" << endl); 00534 Exit(EXIT_DEBUG); 00535 } 00536 #endif 00537 00538 #ifdef CASTOR_VERBOSE 00539 if (m_verbose>=12) 00540 { 00541 string direction = ""; 00542 if (a_direction==FORWARD) direction = "forward"; 00543 else direction = "backward"; 00544 Cout("oProjectionLine::AddVoxelInTOFBin() -> Add voxel index " << a_voxelIndex << " of weight " << a_voxelWeight << 00545 " into TOF bin " << a_TOFBin << " of " << direction << " projector" << endl); 00546 } 00547 #endif 00548 00549 // Different way of doing for each computation strategy 00550 if (m_computationStrategy==IMAGE_COMPUTATION_STRATEGY) 00551 { 00552 // Simply add the contribution to this voxel index 00553 m3p_voxelWeights[a_direction][a_TOFBin][a_voxelIndex] += a_voxelWeight; 00554 } 00555 else if (m_computationStrategy==FIXED_LIST_COMPUTATION_STRATEGY) 00556 { 00557 // No buffer overflow checks for this computation strategy 00558 m3p_voxelIndices[a_direction][a_TOFBin][m2p_currentNbVoxels[a_direction][a_TOFBin]] = a_voxelIndex; 00559 m3p_voxelWeights[a_direction][a_TOFBin][m2p_currentNbVoxels[a_direction][a_TOFBin]] = a_voxelWeight; 00560 m2p_currentNbVoxels[a_direction][a_TOFBin]++; 00561 } 00562 else if (m_computationStrategy==ADAPTATIVE_LIST_COMPUTATION_STRATEGY) 00563 { 00564 // Check if the number of contributing voxels is equal to the allocated size 00565 if (m2p_currentNbVoxels[a_direction][a_TOFBin]==m2p_allocatedNbVoxels[a_direction][a_TOFBin]) 00566 { 00567 // We realloc for one more voxel 00568 m2p_allocatedNbVoxels[a_direction][a_TOFBin]++; 00569 m3p_voxelIndices[a_direction][a_TOFBin] = (INTNB*)realloc(m3p_voxelIndices[a_direction][a_TOFBin],m2p_allocatedNbVoxels[a_direction][a_TOFBin]*sizeof(INTNB)); 00570 m3p_voxelWeights[a_direction][a_TOFBin] = (FLTNB*)realloc(m3p_voxelWeights[a_direction][a_TOFBin],m2p_allocatedNbVoxels[a_direction][a_TOFBin]*sizeof(FLTNB)); 00571 } 00572 // Then register this voxel contribution 00573 m3p_voxelIndices[a_direction][a_TOFBin][m2p_currentNbVoxels[a_direction][a_TOFBin]] = a_voxelIndex; 00574 m3p_voxelWeights[a_direction][a_TOFBin][m2p_currentNbVoxels[a_direction][a_TOFBin]] = a_voxelWeight; 00575 m2p_currentNbVoxels[a_direction][a_TOFBin]++; 00576 } 00577 } 00578 00579 // ===================================================================== 00580 // --------------------------------------------------------------------- 00581 // --------------------------------------------------------------------- 00582 // ===================================================================== 00583 00584 void oProjectionLine::AddVoxel(int a_direction, INTNB a_voxelIndex, FLTNB a_voxelWeight) 00585 { 00586 #ifdef CASTOR_DEBUG 00587 if (!m_initialized) 00588 { 00589 Cerr("***** oProjectionLine::AddVoxel() -> Called while not initialized !" << endl); 00590 Exit(EXIT_DEBUG); 00591 } 00592 #endif 00593 00594 #ifdef CASTOR_VERBOSE 00595 if (m_verbose>=12) 00596 { 00597 string direction = ""; 00598 if (a_direction==FORWARD) direction = "forward"; 00599 else direction = "backward"; 00600 Cout("oProjectionLine::AddVoxel() -> Add voxel index " << a_voxelIndex << " of weight " << a_voxelWeight << 00601 " in " << direction << " projector" << endl); 00602 } 00603 #endif 00604 00605 // No TOF here 00606 int no_tof_bin = 0; 00607 00608 // Different way of doing for each computation strategy 00609 if (m_computationStrategy==IMAGE_COMPUTATION_STRATEGY) 00610 { 00611 // Simply add the contribution to this voxel index 00612 m3p_voxelWeights[a_direction][no_tof_bin][a_voxelIndex] += a_voxelWeight; 00613 } 00614 else if (m_computationStrategy==FIXED_LIST_COMPUTATION_STRATEGY) 00615 { 00616 // No buffer overflow checks for this computation strategy 00617 m3p_voxelIndices[a_direction][no_tof_bin][m2p_currentNbVoxels[a_direction][no_tof_bin]] = a_voxelIndex; 00618 m3p_voxelWeights[a_direction][no_tof_bin][m2p_currentNbVoxels[a_direction][no_tof_bin]] = a_voxelWeight; 00619 m2p_currentNbVoxels[a_direction][no_tof_bin]++; 00620 } 00621 else if (m_computationStrategy==ADAPTATIVE_LIST_COMPUTATION_STRATEGY) 00622 { 00623 // Check if the number of contributing voxels is equal to the allocated size 00624 if (m2p_currentNbVoxels[a_direction][no_tof_bin]==m2p_allocatedNbVoxels[a_direction][no_tof_bin]) 00625 { 00626 // We realloc for one more voxel 00627 m2p_allocatedNbVoxels[a_direction][no_tof_bin]++; 00628 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)); 00629 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)); 00630 } 00631 // Then register this voxel contribution 00632 m3p_voxelIndices[a_direction][no_tof_bin][m2p_currentNbVoxels[a_direction][no_tof_bin]] = a_voxelIndex; 00633 m3p_voxelWeights[a_direction][no_tof_bin][m2p_currentNbVoxels[a_direction][no_tof_bin]] = a_voxelWeight; 00634 m2p_currentNbVoxels[a_direction][no_tof_bin]++; 00635 } 00636 } 00637 00638 // ===================================================================== 00639 // --------------------------------------------------------------------- 00640 // --------------------------------------------------------------------- 00641 // ===================================================================== 00642 00643 FLTNB oProjectionLine::ForwardProject(FLTNB* ap_image) 00644 { 00645 // The forward projection value 00646 FLTNB value = 0.; 00647 // If image is NULL, then assume 1 everywhere 00648 if (ap_image==NULL) 00649 { 00650 // Loop over voxels inside the current TOF bin 00651 for (int v=0; v<m2p_currentNbVoxels[FORWARD][m_currentTOFBin]; v++) 00652 value += m3p_voxelWeights[FORWARD][m_currentTOFBin][v]; 00653 } 00654 // Otherwise, project the image 00655 else 00656 { 00657 if (m_computationStrategy==IMAGE_COMPUTATION_STRATEGY) 00658 { 00659 // Loop over the whole image 00660 for (int v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++) 00661 value += m3p_voxelWeights[FORWARD][m_currentTOFBin][v] * ap_image[v]; 00662 } 00663 else 00664 { 00665 // Loop over voxels inside the current TOF bin 00666 for (int v=0; v<m2p_currentNbVoxels[FORWARD][m_currentTOFBin]; v++) 00667 value += m3p_voxelWeights[FORWARD][m_currentTOFBin][v] 00668 * ap_image[ m3p_voxelIndices[FORWARD][m_currentTOFBin][v] ]; 00669 } 00670 } 00671 // Apply multiplicative correction term 00672 value /= m_multiplicativeCorrection; 00673 // Return the value 00674 return value; 00675 } 00676 00677 // ===================================================================== 00678 // --------------------------------------------------------------------- 00679 // --------------------------------------------------------------------- 00680 // ===================================================================== 00681 00682 FLTNB oProjectionLine::ForwardProjectWithSPECTAttenuation(FLTNB* ap_attenuation, FLTNB* ap_image) 00683 { 00684 #ifdef CASTOR_DEBUG 00685 // This function cannot be used by definition when using the image computation strategy 00686 if (m_computationStrategy==IMAGE_COMPUTATION_STRATEGY) 00687 { 00688 Cerr("***** oProjectionLine::ForwardProjectWithSPECTAttenuation() -> Cannot be used with an image computation strategy over the projection line !" << endl); 00689 Exit(EXIT_FAILURE); // Have to exit 00690 } 00691 // This function cannot be used with an incompatible projector 00692 if (mp_ForwardProjector!=NULL && !mp_ForwardProjector->GetCompatibilityWithSPECTAttenuationCorrection()) 00693 { 00694 Cerr("***** oProjectionLine::ForwardProjectWithSPECTAttenuation() -> Cannot be used with an incompatible projector !" << endl); 00695 Exit(EXIT_FAILURE); // Have to exit 00696 } 00697 // Note that normally, the incompatibilities are checked before launching the reconstruction 00698 #endif 00699 00700 // The forward projection value 00701 FLTNB value = 0.; 00702 00703 // If image is NULL, then assume 1 everywhere in the emission object 00704 if (ap_image==NULL) 00705 { 00706 // Loop over the element in the line 00707 for (int v=0; v<m2p_currentNbVoxels[FORWARD][m_currentTOFBin]; v++) 00708 { 00709 FLTNB atn_sum = 0.0; 00710 // Compute the attenuation of an element belonging to the line 00711 for (int w=v; w<m2p_currentNbVoxels[FORWARD][m_currentTOFBin]; w++) 00712 atn_sum += ap_attenuation[ m3p_voxelIndices[FORWARD][m_currentTOFBin][w] ] * m3p_voxelWeights[FORWARD][m_currentTOFBin][w]; 00713 // Take the inverse exponential to save a division after that (and convert from cm-1 to mm-1) 00714 atn_sum = std::exp( -atn_sum *0.1 ); 00715 value += m3p_voxelWeights[FORWARD][m_currentTOFBin][v] * atn_sum; 00716 } 00717 } 00718 // Otherwise, project the image 00719 else 00720 { 00721 // Loop over the element in the line 00722 for (int v=0; v<m2p_currentNbVoxels[FORWARD][m_currentTOFBin]; v++) 00723 { 00724 FLTNB atn_sum = 0.0; 00725 // Compute the attenuation of an element belonging to the line 00726 for (int w=v; w<m2p_currentNbVoxels[FORWARD][m_currentTOFBin]; w++) 00727 atn_sum += ap_attenuation[ m3p_voxelIndices[FORWARD][m_currentTOFBin][w] ] * m3p_voxelWeights[FORWARD][m_currentTOFBin][w]; 00728 // Take the inverse exponential to save a division after that (and convert from cm-1 to mm-1) 00729 atn_sum = std::exp( -atn_sum *0.1 ); 00730 // Update projection 00731 value += m3p_voxelWeights[FORWARD][m_currentTOFBin][v] * ap_image[ m3p_voxelIndices[FORWARD][m_currentTOFBin][v] ] * atn_sum; 00732 } 00733 } 00734 // Apply multiplicative correction term 00735 value /= m_multiplicativeCorrection; 00736 // Return the value 00737 return value; 00738 } 00739 00740 // ===================================================================== 00741 // --------------------------------------------------------------------- 00742 // --------------------------------------------------------------------- 00743 // ===================================================================== 00744 00745 void oProjectionLine::BackwardProject(FLTNB* ap_image, FLTNB a_value) 00746 { 00747 // Apply multiplicative correction term 00748 a_value /= m_multiplicativeCorrection; 00749 // Image-based strategy 00750 if (m_computationStrategy==IMAGE_COMPUTATION_STRATEGY) 00751 { 00752 // Loop over the whole image 00753 for (int v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++) 00754 ap_image[v] += a_value * m3p_voxelWeights[BACKWARD][m_currentTOFBin][v]; 00755 } 00756 // List-based strategies 00757 else 00758 { 00759 // Loop over voxels inside the current TOF bin 00760 for (int v=0; v<m2p_currentNbVoxels[BACKWARD][m_currentTOFBin]; v++) 00761 ap_image[ m3p_voxelIndices[BACKWARD][m_currentTOFBin][v] ] += a_value * m3p_voxelWeights[BACKWARD][m_currentTOFBin][v]; 00762 } 00763 } 00764 00765 // ===================================================================== 00766 // --------------------------------------------------------------------- 00767 // --------------------------------------------------------------------- 00768 // ===================================================================== 00769 00770 void oProjectionLine::BackwardProjectWithSPECTAttenuation(FLTNB* ap_attenuation, FLTNB* ap_image, FLTNB a_value) 00771 { 00772 #ifdef CASTOR_DEBUG 00773 // This function cannot be used by definition when using the image computation strategy 00774 if (m_computationStrategy==IMAGE_COMPUTATION_STRATEGY) 00775 { 00776 Cerr("***** oProjectionLine::BackwardProjectWithSPECTAttenuation() -> Cannot be used with an image computation strategy over the projection line !" << endl); 00777 Exit(EXIT_FAILURE); // Have to exit 00778 } 00779 // This function cannot be used with an incompatible projector 00780 if (mp_BackwardProjector!=NULL && !mp_BackwardProjector->GetCompatibilityWithSPECTAttenuationCorrection()) 00781 { 00782 Cerr("***** oProjectionLine::BackwardProjectWithSPECTAttenuation() -> Cannot be used with an incompatible projector !" << endl); 00783 Exit(EXIT_FAILURE); // Have to exit 00784 } 00785 // Note that normally, the incompatibilities are checked before launching the reconstruction 00786 #endif 00787 00788 // Apply multiplicative correction term 00789 a_value /= m_multiplicativeCorrection; 00790 00791 // Loop over the element in the line 00792 for (int v=0; v<m2p_currentNbVoxels[BACKWARD][m_currentTOFBin]; v++) 00793 { 00794 FLTNB atn_sum = 0.0; 00795 // Compute the attenuation of an element belonging to the line 00796 for (int w=v; w<m2p_currentNbVoxels[BACKWARD][m_currentTOFBin]; w++) 00797 atn_sum += ap_attenuation[ m3p_voxelIndices[BACKWARD][m_currentTOFBin][w] ] * m3p_voxelWeights[BACKWARD][m_currentTOFBin][w]; 00798 // Take the inverse exponential to save a division after that (and convert from cm-1 to mm-1) 00799 atn_sum = std::exp( -atn_sum *0.1 ); 00800 // Update image 00801 ap_image[ m3p_voxelIndices[BACKWARD][m_currentTOFBin][v] ] += m3p_voxelWeights[BACKWARD][m_currentTOFBin][v] * a_value * atn_sum; 00802 } 00803 } 00804 00805 // ===================================================================== 00806 // --------------------------------------------------------------------- 00807 // --------------------------------------------------------------------- 00808 // ===================================================================== 00809 00810 FLTNB oProjectionLine::ComputeLineIntegral(int a_direction) 00811 { 00812 // Will be the cumulative integral 00813 FLTNB integral = 0.; 00814 // Simply loop over all voxels weights and sum them 00815 for (int v=0; v<m2p_currentNbVoxels[a_direction][m_currentTOFBin]; v++) 00816 integral += m3p_voxelWeights[a_direction][m_currentTOFBin][v]; 00817 // Return the result 00818 return integral; 00819 } 00820