![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class oImageSpace 00004 00005 - separators: X 00006 - doxygen: X 00007 - default initialization: X 00008 - CASTOR_DEBUG: X 00009 - CASTOR_VERBOSE: none 00010 */ 00011 00012 00019 #include "oImageSpace.hh" 00020 #include "vOptimizer.hh" 00021 00022 00023 // ===================================================================== 00024 // --------------------------------------------------------------------- 00025 // --------------------------------------------------------------------- 00026 // ===================================================================== 00027 /* 00028 \brief oImageSpace constructor. 00029 Initialize the member variables to their default values. 00030 */ 00031 oImageSpace::oImageSpace() 00032 { 00033 m_loadedSensitivity = false; 00034 m_loadedAnatomical = false; 00035 m_loadedMask = false; 00036 mp_ID = NULL; 00037 m_verbose = -1; 00038 m_nbBwdImages = 1; 00039 // Set all images pointers to NULL 00040 m4p_image = NULL; 00041 m4p_forwardImage = NULL; 00042 m6p_backwardImage = NULL; 00043 m5p_sensitivity = NULL; 00044 m4p_attenuation = NULL; 00045 m5p_defTmpBackwardImage = NULL; 00046 m4p_defTmpSensitivityImage = NULL; 00047 m4p_outputImage = NULL; 00048 mp_visitedVoxelsImage = NULL; 00049 m2p_projectionImage = NULL; 00050 m4p_anatomicalImage = NULL; 00051 mp_maskImage = NULL; 00052 } 00053 00054 00055 00056 // ===================================================================== 00057 // --------------------------------------------------------------------- 00058 // --------------------------------------------------------------------- 00059 // ===================================================================== 00060 /* 00061 \brief oImageSpace destructor. 00062 */ 00063 oImageSpace::~oImageSpace() {;} 00064 00065 00066 00067 00068 // ===================================================================== 00069 // --------------------------------------------------------------------- 00070 // --------------------------------------------------------------------- 00071 // ===================================================================== 00072 /* 00073 \fn InstantiateImage() 00074 \brief Allocate memory for the main image matrices 00075 */ 00076 void oImageSpace::InstantiateImage() 00077 { 00078 if(m_verbose>=3) Cout("oImageSpace::InstantiateImage ..."<< endl); 00079 00080 // For each temporal dimensions, the number of related images is recovered from the number of intrinsic time basis functions stored in mp_ID 00081 // If no intrinsic temporal basis functions are initialized (standard reconstruction), this is equivalent to the standard calls to GetNbTimeFrames() / GetNbRespGates() / GetNbCardGates() 00082 00083 // First pointer is for the number of time basis functions 00084 m4p_image = new FLTNB***[mp_ID->GetNbTimeBasisFunctions()]; 00085 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 00086 { 00087 // Second pointer is for the number of respiratory basis functions 00088 m4p_image[tbf] = new FLTNB**[mp_ID->GetNbRespBasisFunctions()]; 00089 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 00090 { 00091 // Third pointer is for the number of cardiac basis functions 00092 m4p_image[tbf][rbf] = new FLTNB*[mp_ID->GetNbCardBasisFunctions()]; 00093 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 00094 { 00095 // Fourth pointer is for the 3D space 00096 m4p_image[tbf][rbf][cbf] = new FLTNB[mp_ID->GetNbVoxXYZ()]; 00097 } 00098 } 00099 } 00100 } 00101 00102 00103 00104 00105 00106 // ===================================================================== 00107 // --------------------------------------------------------------------- 00108 // --------------------------------------------------------------------- 00109 // ===================================================================== 00110 /* 00111 \fn DeallocateImage() 00112 \brief Free memory for the main image matrices 00113 */ 00114 void oImageSpace::DeallocateImage() 00115 { 00116 if(m_verbose>=3) Cout("oImageSpace::DeallocateImage ..."<< endl); 00117 00118 // First pointer is for the number of time basis functions 00119 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 00120 { 00121 // Second pointer is for the number of respiratory basis functions 00122 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 00123 { 00124 // Third pointer is for the number of cardiac basis functions 00125 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 00126 { 00127 if (m4p_image[tbf][rbf][cbf]) delete m4p_image[tbf][rbf][cbf]; 00128 } 00129 if (m4p_image[tbf][rbf]) delete[] m4p_image[tbf][rbf]; 00130 } 00131 if (m4p_image[tbf]) delete[] m4p_image[tbf]; 00132 } 00133 if (m4p_image) delete[] m4p_image; 00134 m4p_image = NULL; 00135 } 00136 00137 00138 00139 // ===================================================================== 00140 // --------------------------------------------------------------------- 00141 // --------------------------------------------------------------------- 00142 // ===================================================================== 00143 /* 00144 \fn InstantiateForwardImage() 00145 \brief Allocate memory for the forward image matrices 00146 */ 00147 void oImageSpace::InstantiateForwardImage() 00148 { 00149 if(m_verbose>=3) Cout("oImageSpace::InstantiateForwardImage ..."<< endl); 00150 00151 // For each temporal dimensions, the number of related images is recovered from the number of intrinsic time basis functions stored in mp_ID 00152 // If no intrinsic temporal basis functions are initialized (standard reconstruction), 00153 // this is equivalent to the standard calls to GetNbTimeFrames() / GetNbRespGates() / GetNbCardGates() 00154 00155 // First pointer is for the number of time basis functions 00156 m4p_forwardImage = new FLTNB***[mp_ID->GetNbTimeBasisFunctions()]; 00157 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 00158 { 00159 // Second pointer is for the number of respiratory basis functions 00160 m4p_forwardImage[tbf] = new FLTNB**[mp_ID->GetNbRespBasisFunctions()]; 00161 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 00162 { 00163 // Third pointer is for the number of cardiac basis functions 00164 m4p_forwardImage[tbf][rbf] = new FLTNB*[mp_ID->GetNbCardBasisFunctions()]; 00165 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 00166 { 00167 // Fourth pointer is for the 3D space 00168 m4p_forwardImage[tbf][rbf][cbf] = new FLTNB[mp_ID->GetNbVoxXYZ()]; 00169 } 00170 } 00171 } 00172 } 00173 00174 00175 00176 // ===================================================================== 00177 // --------------------------------------------------------------------- 00178 // --------------------------------------------------------------------- 00179 // ===================================================================== 00180 /* 00181 \fn DeallocateForwardImage() 00182 \brief Free memory for the forward image matrices 00183 */ 00184 void oImageSpace::DeallocateForwardImage() 00185 { 00186 if(m_verbose>=3) Cout("oImageSpace::DeallocateForwardImage ..."<< endl); 00187 00188 // First pointer is for the number of time basis functions 00189 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 00190 { 00191 // Second pointer is for the number of respiratory basis functions 00192 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 00193 { 00194 // Third pointer is for the number of cardiac basis functions 00195 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 00196 { 00197 if (m4p_forwardImage[tbf][rbf][cbf]) delete m4p_forwardImage[tbf][rbf][cbf]; 00198 } 00199 if (m4p_forwardImage[tbf][rbf]) delete[] m4p_forwardImage[tbf][rbf]; 00200 } 00201 if (m4p_forwardImage[tbf]) delete[] m4p_forwardImage[tbf]; 00202 } 00203 if (m4p_forwardImage) delete[] m4p_forwardImage; 00204 m4p_forwardImage = NULL; 00205 } 00206 00207 00208 00209 // ===================================================================== 00210 // --------------------------------------------------------------------- 00211 // --------------------------------------------------------------------- 00212 // ===================================================================== 00213 /* 00214 \fn InstantiateBackwardImage() 00215 \param a_nbBackwardImages : number of backward images required for the optimization algorithm 00216 \brief Allocate memory for the backward image matrices using the number of backward images for the optimization algorithm passed in parameter 00217 */ 00218 void oImageSpace::InstantiateBackwardImage(int a_nbBwdImages) 00219 { 00220 if(m_verbose>=3) Cout("oImageSpace::InstantiateBackwardImage ..."<< endl); 00221 00222 // For each temporal dimensions, the number of related images is recovered from the number of intrinsic time basis functions stored in mp_ID 00223 // If no intrinsic temporal basis functions are initialized (standard reconstruction), this is equivalent to the standard calls to GetNbTimeFrames() / GetNbRespGates() / GetNbCardGates() 00224 00225 // First pointer is for the number of backward images (in case the optimizer does need more than 1) 00226 m_nbBwdImages = a_nbBwdImages; 00227 m6p_backwardImage = new FLTNB*****[m_nbBwdImages]; 00228 for (int img=0; img<m_nbBwdImages; img++) 00229 { 00230 // Second pointer is for the number of threads, in order to have thread-safe backward projections 00231 m6p_backwardImage[img] = new FLTNB****[mp_ID->GetNbThreadsForProjection()]; 00232 for (int th=0; th<mp_ID->GetNbThreadsForProjection(); th++) 00233 { 00234 // The third pointer is for the number of time basis functions 00235 m6p_backwardImage[img][th] = new FLTNB***[mp_ID->GetNbTimeBasisFunctions()]; 00236 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 00237 { 00238 // The fourth pointer is for the number of respiratory basis functions 00239 m6p_backwardImage[img][th][tbf] = new FLTNB**[mp_ID->GetNbRespBasisFunctions()]; 00240 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 00241 { 00242 // The fifth pointer is for the number of cardiac basis functions 00243 m6p_backwardImage[img][th][tbf][rbf] = new FLTNB*[mp_ID->GetNbCardBasisFunctions()]; 00244 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 00245 { 00246 // The sixth pointer is for the 3D space 00247 m6p_backwardImage[img][th][tbf][rbf][cbf] = new FLTNB[mp_ID->GetNbVoxXYZ()]; 00248 } 00249 } 00250 } 00251 } 00252 } 00253 } 00254 00255 00256 00257 // ===================================================================== 00258 // --------------------------------------------------------------------- 00259 // --------------------------------------------------------------------- 00260 // ===================================================================== 00261 /* 00262 \fn DeallocateBackwardImage() 00263 \brief Free memory for the backward image matrices 00264 */ 00265 void oImageSpace::DeallocateBackwardImage() 00266 { 00267 if(m_verbose>=3) Cout("oImageSpace::DeallocateBackwardImage ..."<< endl); 00268 00269 // First pointer is for the number of backward images (in case the optimizer does need more than 1) 00270 for (int img=0; img<m_nbBwdImages; img++) 00271 { 00272 // Second pointer is for the number of threads, in order to have thread-safe backward projections 00273 for (int th=0; th<mp_ID->GetNbThreadsForProjection(); th++) 00274 { 00275 // The third pointer is for the number of time basis functions 00276 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 00277 { 00278 // The fourth pointer is for the number of respiratory basis functions 00279 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 00280 { 00281 // The fifth pointer is for the number of cardiac basis functions 00282 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 00283 { 00284 if (m6p_backwardImage[img][th][tbf][rbf][cbf]) delete m6p_backwardImage[img][th][tbf][rbf][cbf]; 00285 } 00286 if (m6p_backwardImage[img][th][tbf][rbf]) delete m6p_backwardImage[img][th][tbf][rbf]; 00287 } 00288 if (m6p_backwardImage[img][th][tbf]) delete[] m6p_backwardImage[img][th][tbf]; 00289 } 00290 if (m6p_backwardImage[img][th]) delete[] m6p_backwardImage[img][th]; 00291 } 00292 if (m6p_backwardImage[img]) delete[] m6p_backwardImage[img]; 00293 } 00294 if (m6p_backwardImage) delete[] m6p_backwardImage; 00295 m6p_backwardImage = NULL; 00296 } 00297 00298 00299 00300 00301 // ===================================================================== 00302 // --------------------------------------------------------------------- 00303 // --------------------------------------------------------------------- 00304 // ===================================================================== 00305 /* 00306 \fn InstantiateSensitivityImage() 00307 \param a_pathToSensitivityImage : path to the sensitivity image (should be provided only in list-mode reconstruction) 00308 \brief Memory allocation for the sensitivity image matrices 00309 \details Sensitivity image allocation depends on the reconstruction mode : 00310 - list-mode: Sensitivity image has been computed before reconstruction, it is loaded from the path provided in parameter 00311 - histogram: Sensitivity image is calculated on the fly during reconstruction. 00312 First dimension (thread) is only used in histogram mode, as the on-the-fly sensitivity image computation must be thread safe 00313 \return 0 if success, positive value otherwise 00314 */ 00315 void oImageSpace::InstantiateSensitivityImage(const string& a_pathToSensitivityImage) 00316 { 00317 if(m_verbose>=3) Cout("oImageSpace::InstantiateSensitivityImage ..."<< endl); 00318 00319 // ===================================================================================================== 00320 // Case 1: a path to a sensitivity image is provided, meaning that we are in listmode and do not compute 00321 // the sensitivity on-the-fly. In that case, we do not take the multi-threading into account. 00322 // ===================================================================================================== 00323 00324 if (!a_pathToSensitivityImage.empty()) 00325 { 00326 // Allocate for all threads first 00327 m5p_sensitivity = new FLTNB****[mp_ID->GetNbThreadsForProjection()]; 00328 // Allocate the rest only for the first thread 00329 m5p_sensitivity[0] = new FLTNB***[mp_ID->GetNbTimeFrames()]; 00330 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 00331 { 00332 m5p_sensitivity[0][fr] = new FLTNB**[mp_ID->GetNbRespGates()]; 00333 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 00334 { 00335 m5p_sensitivity[0][fr][rg] = new FLTNB*[mp_ID->GetNbCardGates()]; 00336 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 00337 m5p_sensitivity[0][fr][rg][cg] = new FLTNB[mp_ID->GetNbVoxXYZ()]; 00338 } 00339 } 00340 // Make all thread pointers pointing to the first 00341 for (int th=1; th<mp_ID->GetNbThreadsForProjection(); th++) 00342 m5p_sensitivity[th] = m5p_sensitivity[0]; 00343 00344 m_loadedSensitivity = true; 00345 } 00346 00347 // ===================================================================================================== 00348 // Case 2: no sensitivity provided, thus we are in histogram mode and will commpute it on-the-fly. We 00349 // thus need it to be thread-safe, so we allocate for all threads. 00350 // ===================================================================================================== 00351 00352 else 00353 { 00354 // Allocate for all threads first 00355 m5p_sensitivity = new FLTNB****[mp_ID->GetNbThreadsForProjection()]; 00356 // Allocate the rest only for all threads 00357 for (int th=0; th<mp_ID->GetNbThreadsForProjection(); th++) 00358 { 00359 m5p_sensitivity[th] = new FLTNB***[mp_ID->GetNbTimeFrames()]; 00360 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 00361 { 00362 m5p_sensitivity[th][fr] = new FLTNB**[mp_ID->GetNbRespGates()]; 00363 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 00364 { 00365 m5p_sensitivity[th][fr][rg] = new FLTNB*[mp_ID->GetNbCardGates()]; 00366 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 00367 { 00368 m5p_sensitivity[th][fr][rg][cg] = new FLTNB[mp_ID->GetNbVoxXYZ()]; 00369 /* 00370 for (int v=0; v<mp_ID->GetNbVoxXYZ(); v++) 00371 m5p_sensitivity[th][fr][rg][cg][v] = -1.; 00372 * */ 00373 } 00374 } 00375 } 00376 } 00377 m_loadedSensitivity = false; 00378 } 00379 } 00380 00381 00382 00383 00384 // ===================================================================== 00385 // --------------------------------------------------------------------- 00386 // --------------------------------------------------------------------- 00387 // ===================================================================== 00388 /* 00389 \fn InstantiateSensitivityImage() 00390 \param a_pathToSensitivityImage : path to the sensitivity image (should be provided only in list-mode reconstruction) 00391 \brief Free memory for the sensitivity image matrices 00392 \details Sensitivity image deallocation depends on the reconstruction mode (multithreaded in histogram but not in list-mode) 00393 */ 00394 void oImageSpace::DeallocateSensitivityImage() 00395 { 00396 if(m_verbose>=3) Cout("oImageSpace::DeallocateSensitivityImage ..."<< endl); 00397 00398 // If loaded sensitivity (listmode) then fully deallocate only the first pointer 00399 if (m_loadedSensitivity) 00400 { 00401 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 00402 { 00403 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 00404 { 00405 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 00406 { 00407 if (m5p_sensitivity[0][fr][rg][cg]) delete m5p_sensitivity[0][fr][rg][cg]; 00408 } 00409 if (m5p_sensitivity[0][fr][rg]) delete m5p_sensitivity[0][fr][rg]; 00410 } 00411 if (m5p_sensitivity[0][fr]) delete m5p_sensitivity[0][fr]; 00412 } 00413 if (m5p_sensitivity[0]) delete[] m5p_sensitivity[0]; 00414 } 00415 // Otherwise, loop on all threads (histogram) 00416 else 00417 { 00418 for (int th=0; th<mp_ID->GetNbThreadsForProjection(); th++) 00419 { 00420 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 00421 { 00422 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 00423 { 00424 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 00425 { 00426 if (m5p_sensitivity[th][fr][rg][cg]) delete m5p_sensitivity[th][fr][rg][cg]; 00427 } 00428 if (m5p_sensitivity[th][fr][rg]) delete m5p_sensitivity[th][fr][rg]; 00429 } 00430 if (m5p_sensitivity[th][fr]) delete m5p_sensitivity[th][fr]; 00431 } 00432 if (m5p_sensitivity[th]) delete[] m5p_sensitivity[th]; 00433 } 00434 } 00435 // Finally delete the first pointer 00436 if (m5p_sensitivity) delete[] m5p_sensitivity; 00437 m5p_sensitivity = NULL; 00438 } 00439 00440 00441 00442 00443 // ===================================================================== 00444 // --------------------------------------------------------------------- 00445 // --------------------------------------------------------------------- 00446 // ===================================================================== 00447 /* 00448 \fn InitAnatomicalImage() 00449 \param a_pathToImage : path to the anatomical image 00450 \brief Memory allocation and Initialization for the anatomical image matrices 00451 \details Nothing is performed if the path provided in parameter is empty 00452 \return 0 if success, positive value otherwise 00453 */ 00454 int oImageSpace::InitAnatomicalImage(const string& a_pathToImage) 00455 { 00456 if(m_verbose>=3) Cout("oImageSpace::InitAnatomicalImage ..."<< endl); 00457 00458 int nb_atn_frames = mp_ID->GetNbTimeFrames(); 00459 00460 if (!a_pathToImage.empty()) 00461 { 00462 // allocate memory 00463 if (!m_loadedAnatomical) 00464 { 00465 m4p_anatomicalImage = new FLTNB***[nb_atn_frames]; 00466 for (int fr=0; fr<nb_atn_frames; fr++) 00467 { 00468 m4p_anatomicalImage[fr] = new FLTNB**[mp_ID->GetNbRespGates()]; 00469 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 00470 { 00471 m4p_anatomicalImage[fr][rg] = new FLTNB*[mp_ID->GetNbCardGates()]; 00472 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 00473 { 00474 m4p_anatomicalImage[fr][rg][cg] = new FLTNB[mp_ID->GetNbVoxXYZ()]; 00475 } 00476 } 00477 } 00478 } 00479 00480 ifstream image_file(a_pathToImage.c_str(), ios::in|ios::binary); // open file 00481 00482 if(!image_file.is_open()) 00483 { 00484 Cerr("***** oImageSpace::InitAnatomicalImage()-> Error reading file !" << endl); 00485 return 1; 00486 } 00487 else 00488 { 00489 // Interfile image reading (INTF_LERP_ENABLED = interpolation allowed) 00490 if(IntfReadImage(a_pathToImage, m4p_anatomicalImage, mp_ID, m_verbose, INTF_LERP_ENABLED) ) 00491 { 00492 Cerr("***** oImageSpace::InitAnatomicalImage()-> Error reading Interfile : " << a_pathToImage << " !" << endl); 00493 // read failure implies that the anatomical image will never be used, so free all the allocated memory 00494 DeallocateAnatomicalImage(); 00495 return 1; 00496 } 00497 else 00498 { 00499 // flag as loaded 00500 m_loadedAnatomical = true; 00501 } 00502 } 00503 // close file 00504 image_file.close(); 00505 } 00506 00507 // End 00508 return 0; 00509 } 00510 00511 00512 00513 00514 // ===================================================================== 00515 // --------------------------------------------------------------------- 00516 // --------------------------------------------------------------------- 00517 // ===================================================================== 00518 /* 00519 \fn DeallocateAnatomicalImage() 00520 \brief Free memory for the anatomical image matrices 00521 */ 00522 void oImageSpace::DeallocateAnatomicalImage() 00523 { 00524 if(m_verbose>=3) Cout("oImageSpace::DeallocateAnatomicalImage ..."<< endl); 00525 00526 if (m4p_anatomicalImage!=NULL) 00527 { 00528 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 00529 { 00530 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 00531 { 00532 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 00533 { 00534 if (m4p_anatomicalImage[fr][rg][cg]) delete m4p_anatomicalImage[fr][rg][cg]; 00535 } 00536 if (m4p_anatomicalImage[fr][rg]) delete m4p_anatomicalImage[fr][rg]; 00537 } 00538 if (m4p_anatomicalImage[fr]) delete m4p_anatomicalImage[fr]; 00539 } 00540 if (m4p_anatomicalImage) delete[] m4p_anatomicalImage; 00541 m4p_anatomicalImage = NULL; 00542 } 00543 } 00544 00545 00546 00547 00548 // ===================================================================== 00549 // --------------------------------------------------------------------- 00550 // --------------------------------------------------------------------- 00551 // ===================================================================== 00552 /* 00553 \fn InitMaskImage() 00554 \param a_pathToImage : path to the mask image 00555 \brief Memory allocation and Initialization for the mask image matrices 00556 \details Nothing is performed if the path provided in parameter is empty 00557 \return 0 if success, positive value otherwise 00558 */ 00559 int oImageSpace::InitMaskImage(const string& a_pathToImage) 00560 { 00561 if(m_verbose>=3) Cout("oImageSpace::InitMaskImage ..."<< endl); 00562 00563 if (!a_pathToImage.empty()) 00564 { 00565 // allocate memory 00566 if(!m_loadedMask) 00567 { 00568 mp_maskImage = new FLTNB[mp_ID->GetNbVoxXYZ()]; 00569 } 00570 00571 if(m_verbose>=3) Cout("oImageSpace::InitMaskImage ..."<< endl); 00572 00573 ifstream image_file(a_pathToImage.c_str(), ios::in|ios::binary); // open file 00574 00575 if(!image_file.is_open()) 00576 { 00577 Cerr("***** oImageSpace::InitMaskImage()-> Error reading file !" << endl); 00578 return 1; 00579 } 00580 else 00581 { 00582 // Interfile image reading (INTF_LERP_ENABLED = interpolation allowed) 00583 if(IntfReadImage(a_pathToImage, mp_maskImage, mp_ID, m_verbose, INTF_LERP_ENABLED) ) 00584 { 00585 Cerr("***** oImageSpace::InitMaskImage()-> Error reading Interfile : " << a_pathToImage << " !" << endl); 00586 // read failure implies that the mask image will never be used, so free all the allocated memory 00587 DeallocateMaskImage(); 00588 return 1; 00589 } 00590 else 00591 { 00592 // flag as loaded 00593 m_loadedMask = true; 00594 } 00595 } 00596 // close file 00597 image_file.close(); 00598 } 00599 00600 // End 00601 return 0; 00602 } 00603 00604 00605 00606 00607 // ===================================================================== 00608 // --------------------------------------------------------------------- 00609 // --------------------------------------------------------------------- 00610 // ===================================================================== 00611 /* 00612 \fn DeallocateMaskImage() 00613 \brief Free memory for mask image matrix 00614 */ 00615 void oImageSpace::DeallocateMaskImage() 00616 { 00617 if(m_verbose>=3) Cout("oImageSpace::DeallocateMaskImage ..."<< endl); 00618 00619 if (mp_maskImage!=NULL) 00620 { 00621 delete[] mp_maskImage; 00622 mp_maskImage = NULL; 00623 } 00624 00625 } 00626 00627 00628 00629 00630 // ===================================================================== 00631 // --------------------------------------------------------------------- 00632 // --------------------------------------------------------------------- 00633 // ===================================================================== 00634 /* 00635 \fn InstantiateOutputImage() 00636 \brief Instanciate Image matrices dedicated to output writing on disk 00637 \details Additionnal output image matrix is needed in the case the reconstruction uses intrinsic temporal basis functions 00638 In this case, the image matrices are defined in the temporal image basis functions space, therefore requiring an additional step to recover the images in the regular image-space. 00639 If no intrinsic temporal basis functions are used, m4p_outputImage just refers to the address of the image matrix containing the regular image. 00640 */ 00641 void oImageSpace::InstantiateOutputImage() 00642 { 00643 if(m_verbose>=3) Cout("oImageSpace::InstantiateOutputImage ..."<< endl); 00644 00645 // Here, we will allocate the first 3 pointers to the dynamic dimensions 00646 m4p_outputImage = new FLTNB***[mp_ID->GetNbTimeFrames()]; 00647 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 00648 { 00649 m4p_outputImage[fr] = new FLTNB**[mp_ID->GetNbRespGates()]; 00650 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 00651 { 00652 m4p_outputImage[fr][rg] = new FLTNB*[mp_ID->GetNbCardGates()]; 00653 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 00654 { 00655 // For the last pointer (i.e. the voxels), we allocate them only if we strictly need them, which 00656 // means in the case where one of the dynamic dimensions makes internal use of basis functions, 00657 // otherwise, we will just make this output image pointing to the forward image which is useless 00658 // at the moment of computing the output image. 00659 if ( mp_ID->GetTimeStaticFlag() && 00660 mp_ID->GetRespStaticFlag() && 00661 mp_ID->GetCardStaticFlag() ) 00662 { 00663 // In this case, the time frames and basis functions represent the same thing, same for respiratory/cardiac 00664 // gates and their equivalent basis functions 00665 m4p_outputImage[fr][rg][cg] = m4p_forwardImage[fr][rg][cg]; 00666 } 00667 else 00668 { 00669 // We allocate the output image 00670 m4p_outputImage[fr][rg][cg] = new FLTNB[mp_ID->GetNbVoxXYZ()]; 00671 } 00672 } 00673 } 00674 } 00675 } 00676 00677 00678 00679 00680 // ===================================================================== 00681 // --------------------------------------------------------------------- 00682 // --------------------------------------------------------------------- 00683 // ===================================================================== 00684 /* 00685 \fn DeallocateOutputImage() 00686 \brief Free memory for the Image matrices dedicated to output writing on disk 00687 */ 00688 void oImageSpace::DeallocateOutputImage() 00689 { 00690 if(m_verbose>=3) Cout("oImageSpace::DeallocateOutputImage ..."<< endl); 00691 00692 if (m4p_outputImage) 00693 { 00694 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 00695 { 00696 if (m4p_outputImage[fr]) 00697 { 00698 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 00699 { 00700 if (m4p_outputImage[fr][rg]) 00701 { 00702 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 00703 { 00704 // Distinguish between intrinsic use of dynamic basis functions or not (in the latter, 00705 // the m4p_outputImage points to the m4p_forwardImage, so we do not touch it) 00706 if ( (!mp_ID->GetTimeStaticFlag() || 00707 !mp_ID->GetRespStaticFlag() || 00708 !mp_ID->GetCardStaticFlag()) && 00709 m4p_outputImage[fr][rg][cg] ) 00710 { 00711 delete m4p_outputImage[fr][rg][cg]; 00712 } 00713 } 00714 delete[] m4p_outputImage[fr][rg]; 00715 } 00716 } 00717 delete[] m4p_outputImage[fr]; 00718 } 00719 } 00720 delete[] m4p_outputImage; 00721 } 00722 m4p_outputImage = NULL; 00723 } 00724 00725 00726 00727 // ===================================================================== 00728 // --------------------------------------------------------------------- 00729 // --------------------------------------------------------------------- 00730 // ===================================================================== 00731 /* 00732 \fn InstantiateBwdImageForDeformation() 00733 \brief Memory allocation for the buffer backward image required for image-based deformation 00734 */ 00735 void oImageSpace::InstantiateBwdImageForDeformation() 00736 { 00737 if(m_verbose>=3) Cout("oImageSpace::InstantiateBwdImageForDeformation ..."<< endl); 00738 00739 m5p_defTmpBackwardImage = new FLTNB****[m_nbBwdImages]; 00740 00741 for (int img=0; img<m_nbBwdImages; img++) 00742 { 00743 m5p_defTmpBackwardImage[img] = new FLTNB***[mp_ID->GetNbTimeBasisFunctions()]; 00744 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 00745 { 00746 // The fourth pointer is for the number of respiratory basis functions 00747 m5p_defTmpBackwardImage[img][tbf] = new FLTNB**[mp_ID->GetNbRespBasisFunctions()]; 00748 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 00749 { 00750 // The fifth pointer is for the number of cardiac basis functions 00751 m5p_defTmpBackwardImage[img][tbf][rbf] = new FLTNB*[mp_ID->GetNbCardBasisFunctions()]; 00752 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 00753 { 00754 // The sixth pointer is for the 3D space 00755 m5p_defTmpBackwardImage[img][tbf][rbf][cbf] = new FLTNB[mp_ID->GetNbVoxXYZ()]; 00756 } 00757 } 00758 } 00759 } 00760 00761 } 00762 00763 00764 00765 00766 // ===================================================================== 00767 // --------------------------------------------------------------------- 00768 // --------------------------------------------------------------------- 00769 // ===================================================================== 00770 /* 00771 \fn DeallocateBwdImageForDeformation() 00772 \brief Free memory for the buffer backward image required for image-based deformation 00773 */ 00774 void oImageSpace::DeallocateBwdImageForDeformation() 00775 { 00776 if(m_verbose>=3) Cout("oImageSpace::DeallocateBwdImageForDeformation ..."<< endl); 00777 00778 for (int img=0; img<m_nbBwdImages; img++) 00779 { 00780 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 00781 { 00782 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 00783 { 00784 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 00785 { 00786 if (m5p_defTmpBackwardImage[img][tbf][rbf][cbf]) delete m5p_defTmpBackwardImage[img][tbf][rbf][cbf]; 00787 } 00788 if (m5p_defTmpBackwardImage[img][tbf][rbf]) delete m5p_defTmpBackwardImage[img][tbf][rbf]; 00789 } 00790 if (m5p_defTmpBackwardImage[img][tbf]) delete[] m5p_defTmpBackwardImage[img][tbf]; 00791 } 00792 if (m5p_defTmpBackwardImage[img]) delete[] m5p_defTmpBackwardImage[img]; 00793 } 00794 if (m5p_defTmpBackwardImage) delete[] m5p_defTmpBackwardImage; 00795 m5p_defTmpBackwardImage = NULL; 00796 } 00797 00798 00799 00800 00801 00802 // ===================================================================== 00803 // --------------------------------------------------------------------- 00804 // --------------------------------------------------------------------- 00805 // ===================================================================== 00806 /* 00807 \fn InstantiateSensImageForDeformation() 00808 \brief Memory allocation for the buffer sensitivity image required for image-based deformation (only for PET HISTOGRAM mode) 00809 */ 00810 void oImageSpace::InstantiateSensImageForDeformation() 00811 { 00812 if(m_verbose>=3) Cout("oImageSpace::InstantiateSensImageForDeformation ..."<< endl); 00813 00814 m4p_defTmpSensitivityImage = new FLTNB***[mp_ID->GetNbTimeFrames()]; 00815 00816 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 00817 { 00818 m4p_defTmpSensitivityImage[fr] = new FLTNB**[mp_ID->GetNbRespGates()]; 00819 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 00820 { 00821 m4p_defTmpSensitivityImage[fr][rg] = new FLTNB*[mp_ID->GetNbCardGates()]; 00822 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 00823 m4p_defTmpSensitivityImage[fr][rg][cg] = new FLTNB[mp_ID->GetNbVoxXYZ()]; 00824 } 00825 } 00826 } 00827 00828 00829 00830 00831 // ===================================================================== 00832 // --------------------------------------------------------------------- 00833 // --------------------------------------------------------------------- 00834 // ===================================================================== 00835 /* 00836 \fn DeallocateBwdImageForDeformation() 00837 \brief Free memory for the buffer sensitivity image required for image-based deformation 00838 */ 00839 void oImageSpace::DeallocateSensImageForDeformation() 00840 { 00841 if(m_verbose>=3) Cout("oImageSpace::DeallocateSensImageForDeformation ..."<< endl); 00842 00843 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 00844 { 00845 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 00846 { 00847 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 00848 if (m4p_defTmpSensitivityImage[fr][rg][cg]) delete[] m4p_defTmpSensitivityImage[fr][rg][cg]; 00849 00850 if (m4p_defTmpSensitivityImage[fr][rg]) delete[] m4p_defTmpSensitivityImage[fr][rg]; 00851 } 00852 if (m4p_defTmpSensitivityImage[fr]) delete[] m4p_defTmpSensitivityImage[fr]; 00853 } 00854 if (m4p_defTmpSensitivityImage) delete[] m4p_defTmpSensitivityImage; 00855 m4p_defTmpSensitivityImage = NULL; 00856 } 00857 00858 00859 00860 00861 // ===================================================================== 00862 // --------------------------------------------------------------------- 00863 // --------------------------------------------------------------------- 00864 // ===================================================================== 00865 /* 00866 \fn InstantiateVisitedVoxelsImage() 00867 \brief Memory allocation and initialization for the image matrix containing binary information regarding which 3D voxels have been visited during the projection steps. 00868 */ 00869 void oImageSpace::InstantiateVisitedVoxelsImage() 00870 { 00871 if(m_verbose>=3) Cout("oImageSpace::InstantiateVisitedVoxelsImage ..."<< endl); 00872 00873 mp_visitedVoxelsImage = new FLTNB[mp_ID->GetNbVoxXYZ()]; 00874 for (int v=0; v<mp_ID->GetNbVoxXYZ(); v++) 00875 mp_visitedVoxelsImage[v] = 0.; 00876 } 00877 00878 00879 00880 00881 // ===================================================================== 00882 // --------------------------------------------------------------------- 00883 // --------------------------------------------------------------------- 00884 // ===================================================================== 00885 /* 00886 \fn DeallocateVisitedVoxelsImage() 00887 \brief Free memory for the image matrix containing binary information regarding which 3D voxels have been visited during the projection steps. 00888 */ 00889 void oImageSpace::DeallocateVisitedVoxelsImage() 00890 { 00891 if(m_verbose>=3) Cout("oImageSpace::InstantiateVisitedVoxelsImage ..."<< endl); 00892 00893 if (mp_visitedVoxelsImage) delete mp_visitedVoxelsImage; 00894 mp_visitedVoxelsImage = NULL; 00895 } 00896 00897 00898 00899 00900 // ===================================================================== 00901 // --------------------------------------------------------------------- 00902 // --------------------------------------------------------------------- 00903 // ===================================================================== 00904 /* 00905 \fn LMS_DeallocateAttenuationImage() 00906 \brief Free memory for the Attenuation image matrices (for analytical projection or list-mode sensitivity generation) 00907 */ 00908 void oImageSpace::LMS_DeallocateAttenuationImage() 00909 { 00910 if(m_verbose>=3) Cout("oImageSpace::LMS_DeallocateAttenuationImage ..."<< endl); 00911 00912 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 00913 { 00914 for(int rg=0 ; rg<mp_ID->GetSensNbRespGates(); rg++) 00915 { 00916 for(int cg=0 ; cg<mp_ID->GetSensNbCardGates(); cg++) 00917 if(m4p_attenuation[fr][rg][cg] != NULL) delete m4p_attenuation[fr][rg][cg]; 00918 00919 if(m4p_attenuation[fr][rg] != NULL) delete m4p_attenuation[fr][rg]; 00920 } 00921 if(m4p_attenuation[fr] != NULL) delete m4p_attenuation[fr]; 00922 } 00923 if(m4p_attenuation != NULL) delete[] m4p_attenuation; 00924 m4p_attenuation = NULL; 00925 } 00926 00927 00928 00929 00930 // ===================================================================== 00931 // --------------------------------------------------------------------- 00932 // --------------------------------------------------------------------- 00933 // ===================================================================== 00934 /* 00935 \fn InitAttenuationImage() 00936 \param a_pathToAtnImage : path to an existing image 00937 \param a_value : value to initialize each voxel with, if an input image is not provided 00938 \brief Memory allocation and initialisation for the attenuation image using either : - an image provided by the user (a_pathToAtnImage or a_pathToCTImage) 00939 - the default value (=a_value) 00940 \todo If fr/rg/cg > 1 and only one attenuation image is provided, initialize each redundant matrice with pointers to the first one (deal with that here) 00941 \return 0 if success, positive value otherwise 00942 */ 00943 int oImageSpace::InitAttenuationImage(const string& a_pathToAtnImage) 00944 { 00945 if(m_verbose>=3) Cout("oImageSpace::InitAttenuationImage ..."<< endl); 00946 00947 // todo : read the number of dynamic images from Interfile header 00948 // Allocate memory and initialize matrices according to the number of fr/rg/cg 00949 // Provide this info to LMS_LoadAtnImage 00950 00951 // Memory allocation (if not already done) 00952 00953 if(!m4p_attenuation) 00954 { 00955 m4p_attenuation = new FLTNB***[mp_ID->GetNbTimeFrames()]; 00956 00957 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 00958 { 00959 m4p_attenuation[fr] = new FLTNB**[mp_ID->GetSensNbRespGates()]; 00960 00961 for(int rg=0 ; rg<mp_ID->GetSensNbRespGates(); rg++) 00962 { 00963 m4p_attenuation[fr][rg] = new FLTNB*[mp_ID->GetSensNbCardGates()]; 00964 00965 for(int cg=0 ; cg<mp_ID->GetSensNbCardGates(); cg++) 00966 m4p_attenuation[fr][rg][cg] = new FLTNB[mp_ID->GetNbVoxXYZ()]; 00967 } 00968 } 00969 } 00970 00971 00972 // Initialization 00973 00974 if (!a_pathToAtnImage.empty()) // Image initiale 00975 { 00976 if(LoadAttenuationImage(a_pathToAtnImage) ) 00977 { 00978 Cerr("***** oImageSpace::LMS_InitAttenuationImage()-> Error while trying to read file at " << a_pathToAtnImage << endl); 00979 return 1; 00980 } 00981 } 00982 else // Uniform Initialization 00983 { 00984 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 00985 for(int rg=0 ; rg<mp_ID->GetSensNbRespGates() ; rg++) 00986 for(int cg=0 ; cg<mp_ID->GetSensNbCardGates() ; cg++) 00987 for(int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++) 00988 m4p_attenuation[fr][rg][cg][v] = 0.; 00989 } 00990 00991 return 0; 00992 } 00993 00994 00995 00996 00997 // ===================================================================== 00998 // --------------------------------------------------------------------- 00999 // --------------------------------------------------------------------- 01000 // ===================================================================== 01001 /* 01002 \fn LoadAttenuationImage() 01003 \param a_pathToAtnImage : path to an existing image 01004 \brief Load the attenuation image provided by the user in the m2p_attenuation matrix 01005 \todo If fr/rg/cg > 1 and only one attenuation image is provided, initialize each redundant matrice with pointers to the first one (deal with that here) 01006 \return 0 if success, positive value otherwise 01007 */ 01008 int oImageSpace::LoadAttenuationImage(const string& a_pathToImage) 01009 { 01010 if(m_verbose>=3) Cout("oImageSpace::LMS_LoadAttenuationImage ..."<< endl); 01011 01012 ifstream image_file(a_pathToImage.c_str(), ios::in|ios::binary); // Read the corresponding attenuation image 01013 01014 if(!image_file.is_open()) 01015 { 01016 Cerr("***** oImageSpace::LMS_LoadAttenuationImage()-> Error reading file !" << endl); 01017 return 1; 01018 } 01019 else 01020 { 01021 // Interfile image reading (INTF_LERP_ENABLED = interpolation allowed) 01022 if(IntfReadImage(a_pathToImage, m4p_attenuation, mp_ID, m_verbose, INTF_LERP_ENABLED) ) 01023 { 01024 Cerr("***** oImageSpace::PROJ_LoadAttenuationImage()-> Error reading Interfile : " << a_pathToImage << " !" << endl); 01025 return 1; 01026 } 01027 01028 } 01029 01030 image_file.close(); 01031 return 0; 01032 } 01033 01034 01035 01036 01037 01038 01039 // ===================================================================== 01040 // --------------------------------------------------------------------- 01041 // --------------------------------------------------------------------- 01042 // ===================================================================== 01043 /* 01044 \fn InitImage 01045 \param a_pathToInitialImage : path to an existing image 01046 \param a_value : value to initialize each voxel with, if an input image is not provided 01047 \brief Initialize the main image, either using: 01048 - an existing image at path 'a_pathToInitialImage' 01049 - initialize each voxel with 'a_value'. 01050 \return 0 if success, positive value otherwise 01051 */ 01052 int oImageSpace::InitImage(const string& a_pathToInitialImage, FLTNB a_value) 01053 { 01054 if(m_verbose>=3) Cout("oImageSpace::InitImage ..."<< endl); 01055 01056 if (!a_pathToInitialImage.empty()) // Image initiale has been provided 01057 { 01058 if (LoadInitialImage(a_pathToInitialImage) ) 01059 { 01060 Cerr("***** oImageSpace::InitImage()-> Error while trying to read file at " << a_pathToInitialImage << endl); 01061 return 1; 01062 } 01063 } 01064 01065 else // Uniform initialization 01066 { 01067 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 01068 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 01069 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 01070 for (int v=0; v<mp_ID->GetNbVoxXYZ(); v++) 01071 m4p_image[tbf][rbf][cbf][v] = a_value; 01072 } 01073 01074 return 0; 01075 } 01076 01077 01078 01079 01080 // ===================================================================== 01081 // --------------------------------------------------------------------- 01082 // --------------------------------------------------------------------- 01083 // ===================================================================== 01084 /* 01085 \fn LoadInitialImage() 01086 \param a_pathToImage : path to an existing image 01087 \brief Load the initial image provided by the user in the corresponding matrix 01088 \return 0 if success, positive value otherwise 01089 */ 01090 int oImageSpace::LoadInitialImage(const string& a_pathToImage) 01091 { 01092 if(m_verbose>=3) Cout("oImageSpace::LoadInitialImage ..."<< endl); 01093 01094 ifstream image_file; 01095 image_file.open(a_pathToImage.c_str(), ios::in | ios::binary); 01096 01097 if(!image_file.is_open()) 01098 { 01099 Cerr("***** oImageSpace::LoadInitialImage()-> Error reading file!" << endl); 01100 return 1; 01101 } 01102 else 01103 { 01104 // Interfile image reading (INTF_LERP_DISABLED = no interpolation allowed) 01105 // if(IntfReadImage(a_pathToImage, m4p_image, mp_ID, m_verbose, INTF_LERP_DISABLED)) 01106 if(IntfReadImage(a_pathToImage, m4p_image, mp_ID, m_verbose, INTF_LERP_ENABLED)) 01107 { 01108 Cerr("***** oImageSpace::LoadInitialImage()-> Error reading Interfile : " << a_pathToImage << " !" << endl); 01109 return 1; 01110 } 01111 } 01112 image_file.close(); 01113 01114 return 0; 01115 } 01116 01117 01118 01119 01120 // ===================================================================== 01121 // --------------------------------------------------------------------- 01122 // --------------------------------------------------------------------- 01123 // ===================================================================== 01124 /* 01125 \fn InitializationBackwardImage() 01126 \brief Initialize each voxel of the backward images to 0, also for sensitivity if not loaded (estimated on the fly). 01127 */ 01128 void oImageSpace::InitBackwardImage() 01129 { 01130 if(m_verbose>=3) Cout("oImageSpace::InitBackwardImage ..." << endl); 01131 01132 // Reset backward images to 0. 01133 for (int img=0; img<m_nbBwdImages; img++) 01134 { 01135 int th; 01136 #pragma omp parallel for private(th) schedule(static, 1) 01137 for (th=0; th<mp_ID->GetNbThreadsForProjection(); th++) 01138 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 01139 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 01140 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 01141 for (int v=0; v<mp_ID->GetNbVoxXYZ(); v++) 01142 { 01143 m6p_backwardImage[img][th][tbf][rbf][cbf][v] = 0.; 01144 } 01145 } 01146 01147 // If on-the-fly sensitivity, then reset to 0. 01148 if (!m_loadedSensitivity) 01149 { 01150 int th; 01151 #pragma omp parallel for private(th) schedule(static, 1) 01152 for (th=0; th<mp_ID->GetNbThreadsForProjection(); th++) 01153 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 01154 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 01155 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 01156 for (int v=0; v<mp_ID->GetNbVoxXYZ(); v++) 01157 m5p_sensitivity[th][fr][rg][cg][v] = 0.; 01158 } 01159 } 01160 01161 01162 01163 01164 01165 01166 // ===================================================================== 01167 // --------------------------------------------------------------------- 01168 // --------------------------------------------------------------------- 01169 // ===================================================================== 01170 /* 01171 \fn InitSensitivityImage() 01172 \param a_pathToSensitivityImage : path to the sensitivity image (should be provided only in list-mode reconstruction) 01173 \brief Initialization for the sensitivity image matrices 01174 \details Sensitivity image initialization depends on the reconstruction mode : 01175 - list-mode: Sensitivity image has been computed before reconstruction, it is loaded from the path provided in parameter 01176 - histogram: Sensitivity image is calculated on the fly during reconstruction. 01177 First dimension (thread) is only used in histogram mode, as the on-the-fly sensitivity image computation must be thread safe 01178 \return 0 if success, positive value otherwise 01179 */ 01180 int oImageSpace::InitSensitivityImage(const string& a_pathToSensitivityImage) 01181 { 01182 if(m_verbose>=3) Cout("oImageSpace::InitSensitivityImage ..."<< endl); 01183 01184 // ===================================================================================================== 01185 // Case 1: a path to a sensitivity image is provided, meaning that we are in listmode and do not compute 01186 // the sensitivity on-the-fly. In that case, we do not take the multi-threading into account. 01187 // ===================================================================================================== 01188 01189 if (!a_pathToSensitivityImage.empty()) 01190 { 01191 // Open file 01192 ifstream input_file; 01193 input_file.open(a_pathToSensitivityImage.c_str(), ios::binary | ios::in); 01194 // Read file 01195 if (input_file.is_open()) 01196 { 01197 // Interfile image reading (INTF_LERP_DISABLED = no interpolation allowed) 01198 if(IntfReadImage(a_pathToSensitivityImage, m5p_sensitivity[0], mp_ID, m_verbose, INTF_LERP_DISABLED) ) 01199 { 01200 Cerr("***** oImageSpace::InitSensitivityImage()-> Error reading Interfile : " << a_pathToSensitivityImage << " !" << endl); 01201 return 1; 01202 } 01203 input_file.close(); 01204 m_loadedSensitivity = true; 01205 } 01206 else 01207 { 01208 Cerr("***** oImageSpace::InitSensitivityImage() -> Input sensitivity file '" << a_pathToSensitivityImage << "' is missing or corrupted !" << endl); 01209 return 1; 01210 } 01211 } 01212 01213 // ===================================================================================================== 01214 // Case 2: no sensitivity provided, thus we are in histogram mode and will commpute it on-the-fly. We 01215 // thus need it to be thread-safe, so we allocated for all threads. 01216 // ===================================================================================================== 01217 01218 else 01219 { 01220 /* 01221 // Standard initialization 01222 for (int th=0; th<mp_ID->GetNbThreadsForProjection(); th++) 01223 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 01224 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 01225 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 01226 for (int v=0; v<mp_ID->GetNbVoxXYZ(); v++) 01227 m5p_sensitivity[th][fr][rg][cg][v] = -1.; 01228 * */ 01229 // No need to initialize to any value, as it will be done by the InitBackwardImage function which is called at the beginning 01230 // of each subset, and deals with the sensitivity initialization to 0 when m_loadedSensitivity is false 01231 m_loadedSensitivity = false; 01232 } 01233 01234 // End 01235 return 0; 01236 } 01237 01238 01239 01240 01241 01242 // ===================================================================== 01243 // --------------------------------------------------------------------- 01244 // --------------------------------------------------------------------- 01245 // ===================================================================== 01246 /* 01247 \fn InitializationBwdImageForDeformation() 01248 \brief Initialize the buffer backward image dedicated to image-based deformation 01249 */ 01250 void oImageSpace::InitBwdImageForDeformation() 01251 { 01252 if(m_verbose>=3) Cout("oDeformationManager::InitBwdImageForDeformation ..." << endl); 01253 01254 for (int img=0; img<m_nbBwdImages; img++) 01255 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 01256 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 01257 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 01258 for (int v=0; v<mp_ID->GetNbVoxXYZ(); v++) 01259 m5p_defTmpBackwardImage[img][tbf][rbf][cbf][v] = 0.; 01260 } 01261 01262 01263 01264 01265 01266 // ===================================================================== 01267 // --------------------------------------------------------------------- 01268 // --------------------------------------------------------------------- 01269 // ===================================================================== 01270 /* 01271 \fn InitializationSensImageForDeformation() 01272 \brief Initialize the buffer sensitivity image dedicated to image-based deformation, if required (histogram mode, as sensitivity is not loaded) 01273 */ 01274 void oImageSpace::InitSensImageForDeformation() 01275 { 01276 if(m_verbose>=3) Cout("oImageSpace::InitSensImageForDeformation ..."<< endl); 01277 01278 if(m_loadedSensitivity == false) 01279 { 01280 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 01281 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 01282 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 01283 for(int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++) 01284 m4p_defTmpSensitivityImage[fr][rg][cg][v] = 0.; 01285 } 01286 } 01287 01288 01289 01290 01291 01292 // ===================================================================== 01293 // --------------------------------------------------------------------- 01294 // --------------------------------------------------------------------- 01295 // ===================================================================== 01296 /* 01297 \fn InstantiateOutputImage() 01298 \brief Compute output image using the m4p_image matrix and the time/respiratory/cardiac basis functions. Store the result in the m4p_outputImage matrix 01299 \details If time/respiratory/cardiac basis functions have been initialized, this function has no effect. 01300 */ 01301 void oImageSpace::ComputeOutputImage() 01302 { 01303 if(m_verbose>=3) Cout("oImageSpace::ComputeOutputImage ..."<< endl); 01304 01305 // First loop on frames 01306 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 01307 { 01308 // Second loop on respiratory gates 01309 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 01310 { 01311 // Third loop on cardiac gates 01312 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 01313 { 01314 // Reset m4p_outputImage to 0 01315 int v; 01316 #pragma omp parallel for private(v) schedule(guided) 01317 for (v=0; v<mp_ID->GetNbVoxXYZ(); v++) 01318 m4p_outputImage[fr][rg][cg][v] = 0.; 01319 // First loop on time basis functions 01320 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 01321 { 01322 // Get frame/basis coefficient and continue if null 01323 FLTNB time_basis_coef = mp_ID->GetTimeBasisCoefficient(tbf,fr); 01324 if (time_basis_coef==0.) continue; 01325 // Second loop on respiratory basis functions 01326 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 01327 { 01328 // Get resp_gate/basis coefficient and continue if null 01329 FLTNB resp_basis_coef = mp_ID->GetRespBasisCoefficient(rbf,rg); 01330 if (resp_basis_coef==0.) continue; 01331 // Third loop on cardiac basis functions 01332 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 01333 { 01334 // Get card_gate_basis coefficient and continue if null 01335 FLTNB card_basis_coef = mp_ID->GetCardBasisCoefficient(cbf,cg); 01336 if (card_basis_coef==0.) continue; 01337 // Compute global coefficient 01338 FLTNB global_basis_coeff = time_basis_coef * resp_basis_coef * card_basis_coef; 01339 // Loop on voxel with OpenMP (let the chunk size by default as all values are aligned in memory) 01340 #pragma omp parallel for private(v) schedule(guided) 01341 for (v=0; v<mp_ID->GetNbVoxXYZ(); v++) 01342 { 01343 // Add contribution from these basis functions 01344 m4p_outputImage[fr][rg][cg][v] += m4p_image[tbf][rbf][cbf][v] * global_basis_coeff; 01345 } 01346 } 01347 } 01348 } 01349 } 01350 } 01351 } 01352 01353 } 01354 01355 // ===================================================================== 01356 // --------------------------------------------------------------------- 01357 // --------------------------------------------------------------------- 01358 // ===================================================================== 01359 01360 int oImageSpace::ApplyOutputFlip() 01361 { 01362 // If no flip, then return 01363 if ( !mp_ID->GetFlipOutX() && 01364 !mp_ID->GetFlipOutY() && 01365 !mp_ID->GetFlipOutZ() ) return 0; 01366 // Verbose 01367 if (m_verbose>=2) 01368 { 01369 Cout("oImageSpace::ApplyOutputFlip() -> Flip image" << endl); 01370 if (mp_ID->GetFlipOutX()) Cout(" --> Over X" << endl); 01371 if (mp_ID->GetFlipOutY()) Cout(" --> Over Y" << endl); 01372 if (mp_ID->GetFlipOutZ()) Cout(" --> Over Z" << endl); 01373 } 01374 // First loop on frames 01375 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 01376 { 01377 // Second loop on respiratory gates 01378 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 01379 { 01380 // Third loop on cardiac gates 01381 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 01382 { 01383 // Flip over Z 01384 if (mp_ID->GetFlipOutZ()) 01385 { 01386 // Half loop over Z 01387 for (INTNB z_1=0; z_1<mp_ID->GetNbVoxZ()/2; z_1++) 01388 { 01389 // Compute opposite Z 01390 INTNB z_2 = mp_ID->GetNbVoxZ() - 1 - z_1; 01391 // For efficiency 01392 INTNB base_z1 = z_1 * mp_ID->GetNbVoxXY(); 01393 INTNB base_z2 = z_2 * mp_ID->GetNbVoxXY(); 01394 // Loop over Y 01395 for (INTNB y=0; y<mp_ID->GetNbVoxY(); y++) 01396 { 01397 // For efficiency 01398 INTNB base_y = y * mp_ID->GetNbVoxX(); 01399 // Loop over X 01400 for (INTNB x=0; x<mp_ID->GetNbVoxX(); x++) 01401 { 01402 // Compute both indices 01403 INTNB indice_1 = base_z1 + base_y + x; 01404 INTNB indice_2 = base_z2 + base_y + x; 01405 // Switch voxels 01406 FLTNB buffer = m4p_outputImage[fr][rg][cg][indice_1]; 01407 m4p_outputImage[fr][rg][cg][indice_1] = m4p_outputImage[fr][rg][cg][indice_2]; 01408 m4p_outputImage[fr][rg][cg][indice_2] = buffer; 01409 } 01410 } 01411 } 01412 } 01413 // Flip over Y 01414 if (mp_ID->GetFlipOutY()) 01415 { 01416 // Loop over Z 01417 for (INTNB z=0; z<mp_ID->GetNbVoxZ(); z++) 01418 { 01419 // For efficiency 01420 INTNB base_z = z * mp_ID->GetNbVoxXY(); 01421 // Half loop over Y 01422 for (INTNB y_1=0; y_1<mp_ID->GetNbVoxY()/2; y_1++) 01423 { 01424 // Compute opposite Y 01425 INTNB y_2 = mp_ID->GetNbVoxY() - 1 - y_1; 01426 // For efficiency 01427 INTNB base_y1 = y_1 * mp_ID->GetNbVoxX(); 01428 INTNB base_y2 = y_2 * mp_ID->GetNbVoxX(); 01429 // Loop over X 01430 for (INTNB x=0; x<mp_ID->GetNbVoxX(); x++) 01431 { 01432 // Compute both indices 01433 INTNB indice_1 = base_z + base_y1 + x; 01434 INTNB indice_2 = base_z + base_y2 + x; 01435 // Switch voxels 01436 FLTNB buffer = m4p_outputImage[fr][rg][cg][indice_1]; 01437 m4p_outputImage[fr][rg][cg][indice_1] = m4p_outputImage[fr][rg][cg][indice_2]; 01438 m4p_outputImage[fr][rg][cg][indice_2] = buffer; 01439 } 01440 } 01441 } 01442 } 01443 // Flip over X 01444 if (mp_ID->GetFlipOutX()) 01445 { 01446 // Loop over Z 01447 for (INTNB z=0; z<mp_ID->GetNbVoxZ(); z++) 01448 { 01449 // For efficiency 01450 INTNB base_z = z * mp_ID->GetNbVoxXY(); 01451 // Loop over Y 01452 for (INTNB y=0; y<mp_ID->GetNbVoxY(); y++) 01453 { 01454 // For efficiency 01455 INTNB base_y = y * mp_ID->GetNbVoxX(); 01456 // Half loop over X 01457 for (INTNB x_1=0; x_1<mp_ID->GetNbVoxX()/2; x_1++) 01458 { 01459 // Compute opposite X 01460 INTNB x_2 = mp_ID->GetNbVoxX() - 1 - x_1; 01461 // Compute both indices 01462 INTNB indice_1 = base_z + base_y + x_1; 01463 INTNB indice_2 = base_z + base_y + x_2; 01464 // Switch voxels 01465 FLTNB buffer = m4p_outputImage[fr][rg][cg][indice_1]; 01466 m4p_outputImage[fr][rg][cg][indice_1] = m4p_outputImage[fr][rg][cg][indice_2]; 01467 m4p_outputImage[fr][rg][cg][indice_2] = buffer; 01468 } 01469 } 01470 } 01471 } 01472 } 01473 } 01474 } 01475 01476 // Normal end 01477 return 0; 01478 } 01479 01480 // ===================================================================== 01481 // --------------------------------------------------------------------- 01482 // --------------------------------------------------------------------- 01483 // ===================================================================== 01484 01485 int oImageSpace::ApplyOutputFOVMasking() 01486 { 01487 // If the output FOV percent is under 0 (the default value) and number of axial slices to be removed is 0, we do not mask anything 01488 if ( mp_ID->GetFOVOutPercent()<=0. && 01489 mp_ID->GetNbSliceOutMask()==0 ) return 0; 01490 // Verbose 01491 if (m_verbose>=2) 01492 { 01493 Cout("oImageSpace::ApplyOutputFOVMasking() -> Mask output image" << endl); 01494 if (mp_ID->GetFOVOutPercent()>0.) Cout(" --> Mask transaxial FOV outside " << mp_ID->GetFOVOutPercent() << " %" << endl); 01495 if (mp_ID->GetNbSliceOutMask()>0) Cout(" --> Mask " << mp_ID->GetNbSliceOutMask() << " slices from both axial FOV limits" << endl); 01496 } 01497 // ----------------------------------------------- 01498 // Transaxial FOV masking 01499 // ----------------------------------------------- 01500 if (mp_ID->GetFOVOutPercent()>0.) 01501 { 01502 // Precast half the number of voxels over X and Y minus 1 (for efficiency) 01503 FLTNB flt_base_x = 0.5*((FLTNB)(mp_ID->GetNbVoxX()-1)); 01504 FLTNB flt_base_y = 0.5*((FLTNB)(mp_ID->GetNbVoxY()-1)); 01505 // Compute FOV elipse radius over X and Y, then squared 01506 FLTNB squared_radius_x = 0.5 * ((FLTNB)(mp_ID->GetNbVoxX())) * mp_ID->GetVoxSizeX() 01507 * mp_ID->GetFOVOutPercent() / 100.; 01508 squared_radius_x *= squared_radius_x; 01509 FLTNB squared_radius_y = 0.5 * ((FLTNB)(mp_ID->GetNbVoxY())) * mp_ID->GetVoxSizeY() 01510 * mp_ID->GetFOVOutPercent() / 100.; 01511 squared_radius_y *= squared_radius_y; 01512 // We assume that the computation of the distance from the center for a given 01513 // voxel and comparing it with the output FOV percentage costs more than performing 01514 // the loops in an inverse order compared to how the image is stored in memory. 01515 // Thus we begin the loops over X, then Y, then we test and if test passes, we 01516 // do the remaining loop over Z and over all dynamic dimensions. 01517 int x; 01518 #pragma omp parallel for private(x) schedule(guided) 01519 for (x=0; x<mp_ID->GetNbVoxX(); x++) 01520 { 01521 // Compute X distance from image center, then squared 01522 FLTNB squared_distance_x = (((FLTNB)x)-flt_base_x) * mp_ID->GetVoxSizeX(); 01523 squared_distance_x *= squared_distance_x; 01524 // Loop over Y 01525 for (int y=0; y<mp_ID->GetNbVoxY(); y++) 01526 { 01527 // Compute Y distance from image center, then squared 01528 FLTNB squared_distance_y = (((FLTNB)y)-flt_base_y) * mp_ID->GetVoxSizeY(); 01529 squared_distance_y *= squared_distance_y; 01530 // Test if the voxel is inside the FOV elipse, then we skip this voxel 01531 if ( squared_distance_x/squared_radius_x + squared_distance_y/squared_radius_y <= 1. ) continue; 01532 // Loop over Z 01533 for (int z=0; z<mp_ID->GetNbVoxZ(); z++) 01534 { 01535 // Compute global voxel index 01536 INTNB index = z*mp_ID->GetNbVoxXY() + y*mp_ID->GetNbVoxX() + x; 01537 // First loop on frames 01538 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 01539 // Second loop on respiratory gates 01540 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 01541 // Third loop on cardiac gates 01542 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 01543 // Put 0 into this voxel 01544 m4p_outputImage[fr][rg][cg][index] = 0.; 01545 } 01546 } 01547 } 01548 } 01549 // ----------------------------------------------- 01550 // Axial FOV masking 01551 // ----------------------------------------------- 01552 if (mp_ID->GetNbSliceOutMask()>0) 01553 { 01554 INTNB removed_slices = mp_ID->GetNbSliceOutMask(); 01555 // First loop on frames 01556 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 01557 { 01558 // Second loop on respiratory gates 01559 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 01560 { 01561 // Third loop on cardiac gates 01562 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 01563 { 01564 // Mask slices 01565 for (int z=0; z<removed_slices; z++) 01566 { 01567 // First slices 01568 INTNB base_z_first = z*mp_ID->GetNbVoxXY(); 01569 // Loop over Y and X 01570 for (int i=0; i<mp_ID->GetNbVoxXY(); i++) 01571 { 01572 INTNB index = base_z_first + i; 01573 m4p_outputImage[fr][rg][cg][index] = 0.; 01574 } 01575 // Last slices 01576 INTNB base_z_last = (mp_ID->GetNbVoxZ()-1-z)*mp_ID->GetNbVoxXY(); 01577 // Loop over Y and X 01578 for (int i=0; i<mp_ID->GetNbVoxXY(); i++) 01579 { 01580 INTNB index = base_z_last + i; 01581 m4p_outputImage[fr][rg][cg][index] = 0.; 01582 } 01583 } 01584 } 01585 } 01586 } 01587 } 01588 // End 01589 return 0; 01590 } 01591 01592 01593 // ===================================================================== 01594 // --------------------------------------------------------------------- 01595 // --------------------------------------------------------------------- 01596 // ===================================================================== 01597 /* 01598 \fn SaveOutputImage 01599 \param a_iteration : current iteration index 01600 \param a_subset : current number of subsets (or -1 by default) 01601 \brief Call the interfile function to write output image on disk 01602 \return 0 if success, positive value otherwise 01603 */ 01604 int oImageSpace::SaveOutputImage(int a_iteration, int a_subset) 01605 { 01606 if(m_verbose>=3) Cout("oImageSpace::SaveOutputImage ..."<< endl); 01607 01608 // Get the output manager 01609 sOutputManager* p_output_manager = sOutputManager::GetInstance(); 01610 01611 // Interfile writing 01612 // m4p_forwardImage is used here for the special case temporal regularization using intrinsic basis functions is used. 01613 // In this case, m4p_image contains basis function coefficients and not image values. These are computed in ComputeOutputImage() 01614 // In any other case, both matrices contain the image values. 01615 01616 string path_to_img = p_output_manager->GetPathName() + p_output_manager->GetBaseName(); 01617 01618 // Add a suffix for iteration 01619 if (a_iteration >= 0) 01620 { 01621 stringstream ss; ss << a_iteration + 1; 01622 path_to_img.append("_it").append(ss.str()); 01623 } 01624 01625 // Add a suffix for subset (if not negative by default), this means that we save a 'subset' image 01626 if (a_subset >= 0) 01627 { 01628 stringstream ss; ss << a_subset + 1; 01629 path_to_img.append("_ss").append(ss.str()); 01630 } 01631 01632 // We need one "mode" parameter which indicates if we would like to write 1 image file or several 01633 // Adjust functions according to that 01634 if(IntfWriteImgFile(path_to_img, m4p_outputImage, mp_ID, m_verbose) ) 01635 { 01636 Cerr("***** oImageSpace::SaveOutputImage()-> Error writing Interfile of output image !" << endl); 01637 return 1; 01638 } 01639 01640 01641 // Normal end 01642 return 0; 01643 } 01644 01645 01646 01647 01648 01649 // ===================================================================== 01650 // --------------------------------------------------------------------- 01651 // --------------------------------------------------------------------- 01652 // ===================================================================== 01653 /* 01654 \fn SaveDebugImage 01655 \param a_name : output name of the image 01656 \brief Just a debug function dedicated to write any kind of image on disk in raw format, for debugging purposes 01657 */ 01658 void oImageSpace::SaveDebugImage(const string& a_name) 01659 { 01660 #ifdef CASTOR_DEBUG 01661 if(m_verbose>=3) Cout("oImageSpace::SaveDebugImage ..."<< endl); 01662 #endif 01663 01664 ofstream output_file; 01665 output_file.open(a_name.c_str(), ios::binary | ios::out); 01666 01667 //for (int fr=0; fr<ap_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++) 01668 // for (int rg=0; rg<ap_ImageDimensionsAndQuantification->GetNbRespGates(); rg++) 01669 // for (int cg=0; cg<ap_ImageDimensionsAndQuantification->GetNbCardGates(); cg++) 01670 // { 01671 // Write file 01672 output_file.write(reinterpret_cast<char*>(m6p_backwardImage[0][0][0][0][0]), mp_ID->GetNbVoxXYZ()*sizeof(FLTNB)); 01673 // Close file 01674 output_file.close(); 01675 // } 01676 // } 01677 //} 01678 } 01679 01680 01681 01682 01683 // ===================================================================== 01684 // --------------------------------------------------------------------- 01685 // --------------------------------------------------------------------- 01686 // ===================================================================== 01687 /* 01688 \fn PrepareForwardImage 01689 \brief Copy current image matrix in the forward-image buffer matrix 01690 */ 01691 void oImageSpace::PrepareForwardImage() 01692 { 01693 if(m_verbose>=3) Cout("oImageSpace::PrepareForwardImage ..."<< endl); 01694 01695 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 01696 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 01697 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 01698 for (int v=0; v<mp_ID->GetNbVoxXYZ(); v++) 01699 m4p_forwardImage[tbf][rbf][cbf][v] = m4p_image[tbf][rbf][cbf][v]; 01700 } 01701 01702 01703 01704 01705 // ===================================================================== 01706 // --------------------------------------------------------------------- 01707 // --------------------------------------------------------------------- 01708 // ===================================================================== 01709 /* 01710 \fn Reduce 01711 \brief Merge parallel results into the matrix of the backward image matrix of the first thread. Also for MPI. 01712 \todo Choose computation alternative (do not know yet which one is faster...) 01713 */ 01714 void oImageSpace::Reduce() 01715 { 01716 if(m_verbose>=3) Cout("oImageSpace::Reduce ..."<< endl); 01717 01718 #if defined(CASTOR_OMP) || defined(CASTOR_MPI) 01719 if (m_verbose>=3) Cout("oImageSpace::Reduce() -> Merge parallel results" << endl); 01720 #endif 01721 01722 // -------------------------------------------------------------------------------- 01723 // Step 1: merge multi-threads results 01724 // -------------------------------------------------------------------------------- 01725 01726 #ifdef CASTOR_OMP 01727 if (m_verbose>=3) Cout(" --> Over threads ..." << endl); 01728 01729 // Special case here where it appears to be always beneficial to use as many threads 01730 // as the number of threads used for projections. So we set it here, and set it back 01731 // at the end. 01732 omp_set_num_threads(mp_ID->GetNbThreadsForProjection()); 01733 01734 // todo : Choose computation alternative (do not know yet which one is faster...) 01735 int alternative = 2; 01736 01737 // Alternative 1: standard loops based on the pointers hierarchy but mono-thread 01738 if (alternative==1 || mp_ID->GetNbThreadsForImageComputation()==1) 01739 { 01740 for (int img=0; img<m_nbBwdImages; img++) 01741 for (int th=1; th<mp_ID->GetNbThreadsForProjection(); th++) 01742 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 01743 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 01744 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 01745 for (int v=0; v<mp_ID->GetNbVoxXYZ(); v++) 01746 m6p_backwardImage[img][0][tbf][rbf][cbf][v] += m6p_backwardImage[img][th][tbf][rbf][cbf][v]; 01747 } 01748 // Alternative 2: multi-thread merging with the voxels loop at first to be thread safe 01749 // Maybe faster than alternative 1, but the first loop on voxels breaks the memory order... have to try 01750 else if (alternative==2) 01751 { 01752 int v; 01753 #pragma omp parallel for private(v) schedule(guided) 01754 for (v=0; v<mp_ID->GetNbVoxXYZ(); v++) 01755 { 01756 for (int img=0; img<m_nbBwdImages; img++) 01757 for (int th=1; th<mp_ID->GetNbThreadsForProjection(); th++) 01758 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 01759 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 01760 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 01761 m6p_backwardImage[img][0][tbf][rbf][cbf][v] += m6p_backwardImage[img][th][tbf][rbf][cbf][v]; 01762 } 01763 } 01764 // If on-the-fly sensitivity, then do it also for it. 01765 if (!m_loadedSensitivity) 01766 { 01767 if (alternative==1 || mp_ID->GetNbThreadsForImageComputation()==1) 01768 { 01769 for (int th=1; th<mp_ID->GetNbThreadsForProjection(); th++) 01770 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 01771 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 01772 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 01773 for (int v=0; v<mp_ID->GetNbVoxXYZ(); v++) 01774 m5p_sensitivity[0][fr][rg][cg][v] += m5p_sensitivity[th][fr][rg][cg][v]; 01775 } 01776 else if (alternative==2) 01777 { 01778 int v; 01779 #pragma omp parallel for private(v) schedule(guided) 01780 for (v=0; v<mp_ID->GetNbVoxXYZ(); v++) 01781 { 01782 for (int th=1; th<mp_ID->GetNbThreadsForProjection(); th++) 01783 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 01784 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 01785 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 01786 m5p_sensitivity[0][fr][rg][cg][v] += m5p_sensitivity[th][fr][rg][cg][v]; 01787 } 01788 } 01789 } 01790 // Set back the number of threads to the one for image computation 01791 omp_set_num_threads(mp_ID->GetNbThreadsForImageComputation()); 01792 #endif 01793 01794 // -------------------------------------------------------------------------------- 01795 // Step 2: merge multi-instance results (MPI) 01796 // -------------------------------------------------------------------------------- 01797 01798 #ifdef CASTOR_MPI 01799 // We use the MPI_IN_PLACE as the send buffer so that it is the same as the result buffer 01800 if (m_verbose>=3) Cout(" --> Over instances ..." << endl); 01801 01802 // Merge backward images 01803 for (int img=0; img<m_nbBwdImages; img++) 01804 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 01805 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 01806 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 01807 MPI_Allreduce( MPI_IN_PLACE, &m6p_backwardImage[img][0][tbf][rbf][cbf][0], mp_ID->GetNbVoxXYZ(), FLTNBMPI, MPI_SUM, MPI_COMM_WORLD ); 01808 01809 // If on-the-fly sensitivity, then do it also for it 01810 if (!m_loadedSensitivity) 01811 { 01812 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 01813 for (int rg=0; rg<mp_ID->GetNbRespGates(); rg++) 01814 for (int cg=0; cg<mp_ID->GetNbCardGates(); cg++) 01815 MPI_Allreduce( MPI_IN_PLACE, &m5p_sensitivity[0][fr][rg][cg][0], mp_ID->GetNbVoxXYZ(), FLTNBMPI, MPI_SUM, MPI_COMM_WORLD ); 01816 } 01817 #endif 01818 } 01819 01820 01821 01822 01823 01824 // ===================================================================== 01825 // --------------------------------------------------------------------- 01826 // --------------------------------------------------------------------- 01827 // ===================================================================== 01828 /* 01829 \fn CleanNeverVisitedVoxels 01830 \brief Based on the visitedVoxelsImage, clean the never visited voxels in the image. 01831 This function must be called at the end of each iteration 01832 */ 01833 void oImageSpace::CleanNeverVisitedVoxels() 01834 { 01835 if(m_verbose>=3) Cout("oImageSpace::CleanNeverVisitedVoxels ..."<< endl); 01836 01837 int v; 01838 #pragma omp parallel for private(v) schedule(guided) 01839 for (v=0; v<mp_ID->GetNbVoxXYZ(); v++) 01840 { 01841 // Clean this image voxel if the voxel was never visited 01842 if (mp_visitedVoxelsImage[v]==0.) 01843 { 01844 for (int tbf=0; tbf<mp_ID->GetNbTimeBasisFunctions(); tbf++) 01845 for (int rbf=0; rbf<mp_ID->GetNbRespBasisFunctions(); rbf++) 01846 for (int cbf=0; cbf<mp_ID->GetNbCardBasisFunctions(); cbf++) 01847 { 01848 m4p_image[tbf][rbf][cbf][v] = 0.; 01849 } 01850 } 01851 // Reset the image 01852 mp_visitedVoxelsImage[v] = 0.; 01853 } 01854 } 01855 01856 01857 01858 01859 01860 01861 01862 01863 01864 // LIST-MODE SENSITIVITY GENERATION FUNCTIONS 01865 // ===================================================================== 01866 // --------------------------------------------------------------------- 01867 // --------------------------------------------------------------------- 01868 // ===================================================================== 01869 /* 01870 \fn LMS_InstantiateImage() 01871 \brief Allocate memory for the main image matrices (for list-mode sensitivity generation) 01872 \details This function is dedicated to list-mode sensitivity (LMS) generation. 01873 */ 01874 void oImageSpace::LMS_InstantiateImage() 01875 { 01876 if(m_verbose>=3) Cout("oImageSpace::LMS_InstantiateImage ..."<< endl); 01877 01878 m4p_image = new FLTNB***[mp_ID->GetNbTimeFrames()]; 01879 01880 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 01881 { 01882 m4p_image[fr] = new FLTNB**[mp_ID->GetSensNbRespGates()]; 01883 01884 for(int rg=0 ; rg<mp_ID->GetSensNbRespGates() ; rg++) 01885 { 01886 m4p_image[fr][rg] = new FLTNB*[mp_ID->GetSensNbCardGates()]; 01887 01888 for(int cg=0 ; cg<mp_ID->GetSensNbCardGates() ; cg++) 01889 m4p_image[fr][rg][cg] = new FLTNB[mp_ID->GetNbVoxXYZ()]; 01890 } 01891 } 01892 } 01893 01894 01895 01896 01897 // ===================================================================== 01898 // --------------------------------------------------------------------- 01899 // --------------------------------------------------------------------- 01900 // ===================================================================== 01901 /* 01902 \fn LMS_DeallocateImage() 01903 \brief Free memory for the main image matrices (for list-mode sensitivity generation) 01904 \details This function is dedicated to list-mode sensitivity (LMS) generation. 01905 */ 01906 void oImageSpace::LMS_DeallocateImage() 01907 { 01908 if(m_verbose>=3) Cout("oImageSpace::LMS_DeallocateImage ..."<< endl); 01909 01910 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 01911 { 01912 for(int rg=0 ; rg<mp_ID->GetSensNbRespGates() ; rg++) 01913 { 01914 for(int cg=0 ; cg<mp_ID->GetSensNbCardGates() ; cg++) 01915 if(m4p_image[fr][rg][cg] != NULL) delete m4p_image[fr][rg][cg]; 01916 01917 if(m4p_image[fr][rg] != NULL) delete[] m4p_image[fr][rg]; 01918 } 01919 if(m4p_image[fr] != NULL) delete[] m4p_image[fr] ; 01920 } 01921 01922 if(m4p_image != NULL) delete[] m4p_image; 01923 m4p_image = NULL; 01924 } 01925 01926 01927 01928 // ===================================================================== 01929 // --------------------------------------------------------------------- 01930 // --------------------------------------------------------------------- 01931 // ===================================================================== 01932 /* 01933 \fn LMS_InstantiateForwardImage() 01934 \brief Allocate memory for the forward image matrices (for list-mode sensitivity generation) 01935 \details This function is dedicated to list-mode sensitivity (LMS) generation. 01936 */ 01937 void oImageSpace::LMS_InstantiateForwardImage() 01938 { 01939 if(m_verbose>=3) Cout("oImageSpace::LMS_InstantiateForwardImage ..."<< endl); 01940 01941 m4p_forwardImage = new FLTNB***[mp_ID->GetNbTimeFrames()]; 01942 01943 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 01944 { 01945 m4p_forwardImage[fr] = new FLTNB**[mp_ID->GetSensNbRespGates()]; 01946 01947 for(int rg=0 ; rg<mp_ID->GetSensNbRespGates() ; rg++) 01948 { 01949 m4p_forwardImage[fr][rg] = new FLTNB*[mp_ID->GetSensNbCardGates()]; 01950 01951 for(int cg=0 ; cg<mp_ID->GetSensNbCardGates() ; cg++) 01952 { 01953 m4p_forwardImage[fr][rg][cg] = new FLTNB[mp_ID->GetNbVoxXYZ()]; 01954 } 01955 } 01956 01957 } 01958 } 01959 01960 01961 01962 // ===================================================================== 01963 // --------------------------------------------------------------------- 01964 // --------------------------------------------------------------------- 01965 // ===================================================================== 01966 /* 01967 \fn LMS_DeallocateForwardImage() 01968 \brief Free memory for the forward image matrices (for list-mode sensitivity generation) 01969 \details This function is dedicated to list-mode sensitivity (LMS) generation. 01970 */ 01971 void oImageSpace::LMS_DeallocateForwardImage() 01972 { 01973 if(m_verbose>=3) Cout("oImageSpace::LMS_DeallocateForwardImage ..."<< endl); 01974 01975 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 01976 { 01977 for(int rg=0 ; rg<mp_ID->GetSensNbRespGates() ; rg++) 01978 { 01979 for(int cg=0 ; cg<mp_ID->GetSensNbCardGates() ; cg++) 01980 { 01981 if(m4p_forwardImage[fr][rg][cg] != NULL) delete m4p_forwardImage[fr][rg][cg]; 01982 } 01983 01984 if(m4p_forwardImage[fr][rg] != NULL) delete[] m4p_forwardImage[fr][rg]; 01985 } 01986 01987 if(m4p_forwardImage[fr] != NULL) delete[] m4p_forwardImage[fr]; 01988 } 01989 01990 if(m4p_forwardImage != NULL) delete[] m4p_forwardImage; 01991 m4p_forwardImage = NULL; 01992 } 01993 01994 01995 01996 01997 // ===================================================================== 01998 // --------------------------------------------------------------------- 01999 // --------------------------------------------------------------------- 02000 // ===================================================================== 02001 /* 02002 \fn LMS_InstantiateBackwardImage() 02003 \brief Allocate memory for the backward image matrices and initialize them (for list-mode sensitivity generation) 02004 \details This function is dedicated to list-mode sensitivity (LMS) generation. 02005 */ 02006 void oImageSpace::LMS_InstantiateBackwardImage() 02007 { 02008 if(m_verbose>=3) Cout("oImageSpace::LMS_InstantiateBackwardImage ..."<< endl); 02009 02010 m6p_backwardImage = new FLTNB*****[1]; 02011 02012 for(int img=0 ; img<m_nbBwdImages ; img++) 02013 { 02014 m6p_backwardImage[0] = new FLTNB****[mp_ID->GetNbThreadsForProjection()]; 02015 02016 for(int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++) 02017 { 02018 m6p_backwardImage[0][th] = new FLTNB***[mp_ID->GetNbTimeFrames()]; 02019 02020 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 02021 { 02022 m6p_backwardImage[0][th][fr] = new FLTNB**[mp_ID->GetSensNbRespGates()]; 02023 02024 for(int rg=0 ; rg<mp_ID->GetSensNbRespGates() ; rg++) 02025 { 02026 m6p_backwardImage[0][th][fr][rg] = new FLTNB*[mp_ID->GetSensNbCardGates()]; 02027 02028 for(int cg=0 ; cg<mp_ID->GetSensNbCardGates() ; cg++) 02029 { 02030 m6p_backwardImage[0][th][fr][rg][cg] = new FLTNB[mp_ID->GetNbVoxXYZ()]; 02031 02032 for(int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++) 02033 m6p_backwardImage[img][th][fr][rg][cg][v] = 0; 02034 } 02035 } 02036 } 02037 } 02038 } 02039 } 02040 02041 02042 02043 02044 // ===================================================================== 02045 // --------------------------------------------------------------------- 02046 // --------------------------------------------------------------------- 02047 // ===================================================================== 02048 /* 02049 \fn LMS_DeallocateBackwardImage() 02050 \brief Free memory for the backward image matrices (for list-mode sensitivity generation) 02051 \details This function is dedicated to list-mode sensitivity (LMS) generation. 02052 */ 02053 void oImageSpace::LMS_DeallocateBackwardImage() 02054 { 02055 if(m_verbose>=3) Cout("oImageSpace::LMS_DeallocateBackwardImage ..."<< endl); 02056 02057 for(int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++) 02058 { 02059 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 02060 { 02061 for(int rg=0 ; rg<mp_ID->GetSensNbRespGates() ; rg++) 02062 { 02063 for(int cg=0 ; cg<mp_ID->GetSensNbCardGates() ; cg++) 02064 { 02065 if(m6p_backwardImage[0][th][fr][rg][cg] != NULL) delete m6p_backwardImage[0][th][fr][rg][cg]; 02066 } 02067 02068 if(m6p_backwardImage[0][th][fr][rg] != NULL) delete m6p_backwardImage[0][th][fr][rg]; 02069 } 02070 if(m6p_backwardImage[0][th][fr] != NULL) delete[] m6p_backwardImage[0][th][fr]; 02071 } 02072 if(m6p_backwardImage[0][th] != NULL) delete[] m6p_backwardImage[0][th]; 02073 } 02074 02075 if(m6p_backwardImage[0] != NULL) delete[] m6p_backwardImage[0]; 02076 if(m6p_backwardImage != NULL) delete[] m6p_backwardImage; 02077 m6p_backwardImage = NULL; 02078 } 02079 02080 02081 02082 02083 // ===================================================================== 02084 // --------------------------------------------------------------------- 02085 // --------------------------------------------------------------------- 02086 // ===================================================================== 02087 /* 02088 \fn LMS_InstantiateSensitivityImage() 02089 \brief Allocate memory for the sensitivity image matrices (for list-mode sensitivity generation) 02090 \details This function is dedicated to list-mode sensitivity (LMS) generation. 02091 */ 02092 void oImageSpace::LMS_InstantiateSensitivityImage() 02093 { 02094 if(m_verbose>=3) Cout("oImageSpace::LMS_InstantiateSensitivityImage ..."<< endl); 02095 02096 m5p_sensitivity = new FLTNB****[1]; 02097 m5p_sensitivity[0] = new FLTNB***[mp_ID->GetNbTimeFrames()]; 02098 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 02099 { 02100 m5p_sensitivity[0][fr] = new FLTNB**[mp_ID->GetSensNbRespGates()]; 02101 for (int rg=0; rg<mp_ID->GetSensNbRespGates(); rg++) 02102 { 02103 m5p_sensitivity[0][fr][rg] = new FLTNB*[mp_ID->GetSensNbCardGates()]; 02104 for (int cg=0; cg<mp_ID->GetSensNbCardGates(); cg++) 02105 m5p_sensitivity[0][fr][rg][cg] = new FLTNB[mp_ID->GetNbVoxXYZ()]; 02106 } 02107 } 02108 } 02109 02110 02111 02112 02113 // ===================================================================== 02114 // --------------------------------------------------------------------- 02115 // --------------------------------------------------------------------- 02116 // ===================================================================== 02117 /* 02118 \fn LMS_DeallocateSensitivityImage() 02119 \brief Free memory for the sensitivity image matrices (for list-mode sensitivity generation) 02120 \details This function is dedicated to list-mode sensitivity (LMS) generation. 02121 */ 02122 void oImageSpace::LMS_DeallocateSensitivityImage() 02123 { 02124 if(m_verbose>=3) Cout("oImageSpace::LMS_DeallocateSensitivityImage ..."<< endl); 02125 02126 for (int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++) 02127 { 02128 for (int rg=0; rg<mp_ID->GetSensNbRespGates(); rg++) 02129 { 02130 for (int cg=0; cg<mp_ID->GetSensNbCardGates(); cg++) 02131 { 02132 if(m5p_sensitivity[0][fr][rg][cg] != NULL) delete m5p_sensitivity[0][fr][rg][cg]; 02133 } 02134 if(m5p_sensitivity[0][fr][rg] != NULL)delete m5p_sensitivity[0][fr][rg]; 02135 } 02136 if(m5p_sensitivity[0][fr] != NULL)delete m5p_sensitivity[0][fr]; 02137 } 02138 if(m5p_sensitivity[0] != NULL) delete[] m5p_sensitivity[0]; 02139 02140 if(m5p_sensitivity != NULL) delete[] m5p_sensitivity; 02141 m5p_sensitivity = NULL; 02142 } 02143 02144 02145 02146 02147 // ===================================================================== 02148 // --------------------------------------------------------------------- 02149 // --------------------------------------------------------------------- 02150 // ===================================================================== 02151 /* 02152 \fn LMS_CopyAtnToImage() 02153 \brief Copy the attenuation image contained in the 'm2p_attenuation' matrix inside the m2p_image matrix (for list-mode sensitivity generation) 02154 \details This function is dedicated to list-mode sensitivity (LMS) generation. 02155 */ 02156 void oImageSpace::LMS_CopyAtnToImage() 02157 { 02158 if(m_verbose>=3) Cout("oImageSpace::LMS_CopyAtnToImage ..."<< endl); 02159 02160 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 02161 for(int rg=0 ; rg<mp_ID->GetSensNbRespGates() ; rg++) 02162 for(int cg=0 ; cg<mp_ID->GetSensNbCardGates() ; cg++) 02163 for(int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++) 02164 { 02165 m4p_image[fr][rg][cg][v] = m4p_attenuation[fr][rg][cg][v]; 02166 } 02167 02168 } 02169 02170 02171 // ===================================================================== 02172 // --------------------------------------------------------------------- 02173 // --------------------------------------------------------------------- 02174 // ===================================================================== 02175 /* 02176 \fn LMS_CopyAtnToForwardImage() 02177 \brief Copy the attenuation image contained in the 'm2p_attenuation' matrix inside the m2p_forwardImage matrix (for list-mode sensitivity generation) 02178 \details This function is dedicated to list-mode sensitivity (LMS) generation. 02179 */ 02180 void oImageSpace::LMS_CopyAtnToForwardImage() 02181 { 02182 if(m_verbose>=3) Cout("oImageSpace::LMS_CopyAtnToForwardImage ..."<< endl); 02183 02184 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 02185 for(int rg=0 ; rg<mp_ID->GetSensNbRespGates() ; rg++) 02186 for(int cg=0 ; cg<mp_ID->GetSensNbCardGates() ; cg++) 02187 for(int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++) 02188 { 02189 m4p_forwardImage[fr][rg][cg][v] = m4p_attenuation[fr][rg][cg][v]; 02190 } 02191 02192 } 02193 02194 // ===================================================================== 02195 // --------------------------------------------------------------------- 02196 // --------------------------------------------------------------------- 02197 // ===================================================================== 02198 /* 02199 \fn LMS_CopyBackwardToSensitivityImage() 02200 \brief Copy the backward image containing the result of the sensitivity back-projection, into the sensitivity matrix. 02201 \details This function is dedicated to list-mode sensitivity (LMS) generation, it is useful to apply convolution and to save the image. 02202 */ 02203 void oImageSpace::LMS_CopyBackwardToSensitivity() 02204 { 02205 if(m_verbose>=3) Cout("oImageSpace::LMS_CopyBackwardToSensitivity ..."<< endl); 02206 for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 02207 for (int rg=0 ; rg<mp_ID->GetSensNbRespGates() ; rg++) 02208 for (int cg=0 ; cg<mp_ID->GetSensNbCardGates() ; cg++) 02209 for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++) 02210 m5p_sensitivity[0][fr][rg][cg][v] = m6p_backwardImage[0][0][fr][rg][cg][v]; 02211 } 02212 02213 02214 // ===================================================================== 02215 // --------------------------------------------------------------------- 02216 // --------------------------------------------------------------------- 02217 // ===================================================================== 02218 /* 02219 \fn LMS_PrepareForwardImage() 02220 \brief Copy current image in forward-image buffer (for list-mode sensitivity generation) 02221 \details This function is dedicated to list-mode sensitivity (LMS) generation. 02222 */ 02223 void oImageSpace::LMS_PrepareForwardImage() 02224 { 02225 if(m_verbose>=3) Cout("oImageSpace::LMS_PrepareForwardImage ..."<< endl); 02226 02227 for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++) 02228 for (int rg=0 ; rg<mp_ID->GetSensNbRespGates() ; rg++) 02229 for (int cg=0 ; cg<mp_ID->GetSensNbCardGates() ; cg++) 02230 for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++) 02231 m4p_forwardImage[fr][rg][cg][v] = m4p_image[fr][rg][cg][v]; 02232 } 02233 02234 02235 02236 02237 // ===================================================================== 02238 // --------------------------------------------------------------------- 02239 // --------------------------------------------------------------------- 02240 // ===================================================================== 02241 /* 02242 \fn LMS_Reduce() 02243 \param a_fr : time frame 02244 \param a_rg : respiratory gate 02245 \param a_cg : cardiac gate 02246 \brief Merge parallel results into the matrix of the backward image matrix of the first thread for the specific frame / respiratory gate / cardiac gate (for list-mode sensitivity generation) 02247 \details This function is dedicated to list-mode sensitivity (LMS) generation. 02248 */ 02249 void oImageSpace::LMS_Reduce(int a_fr, int a_rg, int a_cg) 02250 { 02251 if(m_verbose>=3) Cout("oImageSpace::LMS_Reduce ..."<< endl); 02252 02253 for (int th=1 ; th<mp_ID->GetNbThreadsForProjection() ; th++) 02254 for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++) 02255 m6p_backwardImage[0][0][a_fr][a_rg][a_cg][v] += m6p_backwardImage[0][th][a_fr][a_rg][a_cg][v]; 02256 02257 } 02258 02259 02260 02261 02262 // ===================================================================== 02263 // --------------------------------------------------------------------- 02264 // --------------------------------------------------------------------- 02265 // ===================================================================== 02266 /* 02267 \fn LMS_SaveSensitivityImage 02268 \param a_pathToSensitivityImage : path to the sensitivity image (should be provided only in list-mode reconstruction) 02269 \param ap_DeformationManager : Pointer to the deformation manager objet (required to retrieve the number of gates in the sensitivity image) 02270 \brief Call the interfile function to write the sensitivity image on disk 02271 \details If image deformation is enabled for respiratory/cardiac gated data, the gated images are summed up into one image and normalize 02272 \return 0 if success, positive value otherwise 02273 \todo Interfile management : if the number of sensitivity image to write is different to the actual number of image 02274 as it could be the case for dynamic imaging, if one sensitivity image is required for the whole dynamic serie, 02275 we will have to give this information to IntfWriteImgFile(), otherwise the Interfile functions will try to write several times the image (leading to segfault or error) 02276 */ 02277 int oImageSpace::LMS_SaveSensitivityImage(const string& a_pathToSensitivityImage, oDeformationManager* ap_DeformationManager) 02278 { 02279 if(m_verbose>=3) Cout("oImageSpace::LMS_SaveSensitivityImage ..."<< endl); 02280 02281 //sOutputManager* p_output_manager = sOutputManager::GetInstance(); 02282 02283 // According to the respiratory motion algorithm, images relatives to resp gates, cardiac gates, or both, have to be summed up and normalize into one image. 02284 // The resulting image will be summed up into the sensitivity matrix 02285 02286 int nb_reco_card_images = ap_DeformationManager->GetNbSensImagesCardDeformation(mp_ID->GetSensNbCardGates()); 02287 02288 if (ap_DeformationManager->GetNbSensImagesCardDeformation(mp_ID->GetSensNbCardGates()) == 0) // Summation of the card images 02289 { 02290 nb_reco_card_images = 1; 02291 FLTNB card_normalization = mp_ID->GetSensNbCardGates(); 02292 02293 for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++) 02294 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames(); fr++) 02295 for (int rg=0 ; rg<mp_ID->GetSensNbRespGates() ; rg++) 02296 { 02297 for (int cg=1 ; cg<mp_ID->GetSensNbCardGates() ; cg++) 02298 m5p_sensitivity[0][fr][rg][0][v] += m5p_sensitivity[0][fr][rg][cg][v]; 02299 m5p_sensitivity[0][fr][rg][0][v] /= card_normalization; 02300 } 02301 } 02302 02303 if (ap_DeformationManager->GetNbSensImagesRespDeformation(mp_ID->GetSensNbRespGates()) == 0) // Summation of the resp images 02304 { 02305 FLTNB resp_normalization = mp_ID->GetSensNbRespGates(); 02306 02307 for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++) 02308 for(int fr=0 ; fr<mp_ID->GetNbTimeFrames(); fr++) 02309 for (int cg=0 ; cg<nb_reco_card_images ; cg++) // Dual gating : Be sure to get the right number of card images 02310 { 02311 for (int rg=1 ; rg<mp_ID->GetSensNbRespGates() ; rg++) 02312 m5p_sensitivity[0][fr][0][cg][v] += m5p_sensitivity[0][fr][rg][cg][v]; 02313 m5p_sensitivity[0][fr][0][cg][v] /= resp_normalization; 02314 } 02315 } 02316 02317 if (IntfWriteImgFile(a_pathToSensitivityImage, m5p_sensitivity[0], mp_ID, m_verbose) ) 02318 { 02319 Cerr("***** oImageSpace::LMS_SaveSensitivityImage()-> Error writing Interfile of output image !" << endl); 02320 return 1; 02321 } 02322 02323 return 0; 02324 } 02325 02326 02327 // PROJECTION SCRIPT FUNCTIONS 02328 // ===================================================================== 02329 // --------------------------------------------------------------------- 02330 // --------------------------------------------------------------------- 02331 // ===================================================================== 02332 /* 02333 \fn PROJ_InstantiateProjectionImage() 02334 \param a_nbProjs : a number of projection slices in the projection 02335 \param a_nbPixels : a total number of pixels in the projection slices 02336 \brief Instanciate and initialize projection image for analytical projection 02337 \details This function is currently only dedicated to SPECT projection, and used by the analytical projection script 02338 \todo support for PET (requires projection format) 02339 */ 02340 void oImageSpace::PROJ_InstantiateProjectionImage(int a_nbProjs, int a_nbPixels) 02341 { 02342 if(m_verbose>=3) Cout("oImageSpace::PROJ_InstantiateProjectionImage ..."<< endl); 02343 02344 m2p_projectionImage = new FLTNB*[a_nbProjs]; 02345 02346 for(int p=0 ; p<a_nbProjs ; p++) 02347 { 02348 m2p_projectionImage[p] = new FLTNB[a_nbPixels]; 02349 for(int px=0 ; px<a_nbPixels ; px++) 02350 m2p_projectionImage[p][px] = 0; 02351 } 02352 } 02353 02354 02355 02356 02357 // ===================================================================== 02358 // --------------------------------------------------------------------- 02359 // --------------------------------------------------------------------- 02360 // ===================================================================== 02361 /* 02362 \fn PROJ_DeallocateProjectionImage() 02363 \param a_nbProjs : a number of projection slices in the projection 02364 \brief Free memory for the projection image for analytical projection 02365 \details This function is currently only dedicated to SPECT projection, and used by the analytical projection script 02366 \todo support for PET (requires projection format) 02367 */ 02368 void oImageSpace::PROJ_DeallocateProjectionImage(int a_nbProjs) 02369 { 02370 if(m_verbose>=3) Cout("oImageSpace::PROJ_DeallocateProjectionImage ..."<< endl); 02371 02372 for(int p=0 ; p<a_nbProjs ; p++) 02373 { 02374 delete[] m2p_projectionImage[p]; 02375 } 02376 02377 delete[] m2p_projectionImage; 02378 m2p_projectionImage = NULL; 02379 } 02380 02381 02382 02383 02384 // ===================================================================== 02385 // --------------------------------------------------------------------- 02386 // --------------------------------------------------------------------- 02387 // ===================================================================== 02388 /* 02389 \fn PROJ_InitImage() 02390 \param a_pathToInitialImage : path to the image to project 02391 \brief Load the initial image for the analytical projection 02392 \return 0 if success, positive value otherwise 02393 */ 02394 int oImageSpace::PROJ_InitImage(const string& a_pathToInitialImage) 02395 { 02396 if(m_verbose>=3) Cout("oImageSpace::PROJ_InitImage ..."<< endl); 02397 02398 if (!a_pathToInitialImage.empty()) // Image initiale 02399 { 02400 if(PROJ_LoadInitialImage(a_pathToInitialImage) ) 02401 { 02402 Cerr("***** oImageSpace::PROJ_InitImage()-> Error while trying to read file at " << a_pathToInitialImage << endl); 02403 return 1; 02404 } 02405 } 02406 else 02407 { 02408 { 02409 Cerr("***** oImageSpace::PROJ_InitImage()-> No projected image provided ! " << endl); 02410 return 1; 02411 } 02412 } 02413 return 0; 02414 } 02415 02416 02417 02418 02419 // ===================================================================== 02420 // --------------------------------------------------------------------- 02421 // --------------------------------------------------------------------- 02422 // ===================================================================== 02423 /* 02424 \fn PROJ_LoadInitialImage() 02425 \param a_pathToImage : path to the image to project 02426 \brief Load the initial image for the analytical projection 02427 \return 0 if success, positive value otherwise 02428 */ 02429 int oImageSpace::PROJ_LoadInitialImage(const string& a_pathToImage) 02430 { 02431 if(m_verbose>=3) Cout("oImageSpace::PROJ_LoadInitialImage ..."<< endl); 02432 02433 ifstream image_file; 02434 image_file.open(a_pathToImage.c_str(), ios::in | ios::binary); 02435 02436 if(!image_file.is_open()) 02437 { 02438 Cerr("***** oImageSpace::PROJ_LoadInitialImage()-> Error reading file !" << endl); 02439 return 1; 02440 } 02441 else 02442 { 02443 // Interfile image reading (INTF_LERP_DISABLED = no interpolation allowed) 02444 if(IntfReadImage(a_pathToImage, m4p_image, mp_ID, m_verbose, INTF_LERP_DISABLED)) 02445 { 02446 Cerr("***** oImageSpace::PROJ_LoadInitialImage()-> Error reading Interfile : " << a_pathToImage << " !" << endl); 02447 return 1; 02448 } 02449 02450 } 02451 image_file.close(); 02452 return 0; 02453 } 02454 02455 02456 02457 02458 02459 // ===================================================================== 02460 // --------------------------------------------------------------------- 02461 // --------------------------------------------------------------------- 02462 // ===================================================================== 02463 /* 02464 \fn PROJ_SaveProjectionImage 02465 \brief Save an image of the projected data 02466 \return 0 if success, positive value otherwise 02467 \todo Support for PET projeted data 02468 implement interfile file recovery in a scanner class function 02469 */ 02470 int oImageSpace::PROJ_SaveProjectionImage() 02471 { 02472 if(m_verbose>=3) Cout("oImageSpace::PROJ_SaveProjectionImage ..."<< endl); 02473 02474 string img_name = "_ProjectionImage"; 02475 02476 // Initialize interfile fields structure 02477 Intf_fields Img_fields; 02478 IntfKeySetFieldsOutput(&Img_fields , mp_ID); 02479 02480 Img_fields.process_status = "acquired"; 02481 02482 // Get specific info about projection 02483 sScannerManager* p_ScanMgr = sScannerManager::GetInstance(); 02484 02485 // Recover the values in the interfile structure. 02486 // todo : This step for PET projection 02487 // todo : Implement this step in scanner class functions 02488 // todo : For spect, deal with systems with several detector heads (currently assume one head) 02489 if(p_ScanMgr->GetScannerType() == 2) 02490 { 02491 uint16_t nb_projs, nb_heads; 02492 uint16_t nb_bins[2]; 02493 FLTNB pix_size[2]; 02494 FLTNB* p_angles; 02495 FLTNB* p_radius; 02496 int dir_rot; 02497 p_ScanMgr->GetSPECTSpecificParameters(&nb_projs, 02498 &nb_heads, 02499 nb_bins, 02500 pix_size, 02501 p_angles, 02502 p_radius, 02503 &dir_rot); 02504 02505 Img_fields.nb_projections = nb_projs; 02506 Img_fields.nb_detector_heads = nb_heads; 02507 Img_fields.mtx_size[0] = nb_bins[0]; 02508 Img_fields.mtx_size[1] = nb_bins[1]; 02509 Img_fields.vox_size[0] = pix_size[0]; 02510 Img_fields.vox_size[1] = pix_size[1]; 02511 Img_fields.extent_rotation = 360; // 360 by default 02512 // Check rotation direction from the two first angles 02513 Img_fields.direction_rotation = (dir_rot == GEO_ROT_CCW)? 02514 "CCW" : 02515 "CW"; 02516 Img_fields.first_angle = p_angles[0]; 02517 02518 // Note: In Interfile v3.3, doesn't seem to have a field to provide 02519 // individually each angle and Center of rotation. 02520 // Looks like it was planned for ulterior version, at least for the 02521 // center of rotation 02522 // For now, we just write each projection angle and radius as an array key 02523 02524 string angles_str = "{"; 02525 string radius_str = "{"; 02526 02527 // Flag to check if the radius is similar for each projection angle 02528 bool has_single_radius = true; 02529 02530 for(uint16_t p=0 ; p<nb_projs ; p++) 02531 { 02532 stringstream ss_a, ss_r; 02533 // projection angles 02534 ss_a << p_angles[p]; 02535 angles_str.append(ss_a.str()); 02536 (p<nb_projs-1) ? angles_str.append(",") : angles_str.append("}"); 02537 02538 // radius 02539 ss_r << p_radius[p]; 02540 radius_str.append(ss_r.str()); 02541 (p<nb_projs-1) ? radius_str.append(",") : radius_str.append("}"); 02542 if(p_radius[p] != p_radius[0]) 02543 has_single_radius = false; 02544 02545 // Some interfile editors struggle with long line, perhaps because 02546 // they are stored in char[256] or something like that 02547 // Hence, break line after a certain number of elements 02548 if((p+1)%30 == 0) 02549 { 02550 angles_str.append("\n"); 02551 radius_str.append("\n"); 02552 } 02553 } 02554 02555 // If there is one common radius, write its value in the related string 02556 if(has_single_radius) 02557 { 02558 stringstream ss; 02559 ss << p_radius[0]; 02560 radius_str = ss.str(); 02561 } 02562 02563 Img_fields.projection_angles = angles_str; 02564 Img_fields.radius = radius_str; 02565 02566 02567 // Common code to all modalities 02568 sOutputManager* p_output_manager = sOutputManager::GetInstance(); 02569 string img_file_name = p_output_manager->GetPathName() + p_output_manager->GetBaseName(); 02570 img_file_name.append("_ProjImage"); 02571 02572 if(IntfWriteProjFile(img_file_name, m2p_projectionImage, mp_ID, Img_fields, m_verbose) ) 02573 { 02574 Cerr("***** oImageSpace::PROJ_SaveProjectionImage()-> Error writing Interfile of output image !" << endl); 02575 return 1; 02576 } 02577 } 02578 return 0; 02579 } 02580 02581