![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class oImageDimensionsAndQuantification 00004 00005 - separators: X 00006 - doxygen: 00007 - default initialization: X 00008 - CASTOR_DEBUG: 00009 - CASTOR_VERBOSE: 00010 */ 00011 00018 #include "oImageDimensionsAndQuantification.hh" 00019 00020 // ===================================================================== 00021 // --------------------------------------------------------------------- 00022 // --------------------------------------------------------------------- 00023 // ===================================================================== 00024 00025 oImageDimensionsAndQuantification::oImageDimensionsAndQuantification() 00026 { 00027 // Set all members to default values 00028 m_nbVoxX = -1; 00029 m_nbVoxY = -1; 00030 m_nbVoxZ = -1; 00031 m_nbVoxXY = -1; 00032 m_nbVoxXYZ = -1; 00033 m_voxSizeX = -1.; 00034 m_voxSizeY = -1.; 00035 m_voxSizeZ = -1.; 00036 m_fovSizeX = -1.; 00037 m_fovSizeY = -1.; 00038 m_fovSizeZ = -1.; 00039 m_fovOutPercent = 0.; 00040 m_flipOutX = false; 00041 m_flipOutY = false; 00042 m_flipOutZ = false; 00043 m_offsetX = 0.; 00044 m_offsetY = 0.; 00045 m_offsetZ = 0.; 00046 m_frameList = ""; 00047 m_nbTimeFrames = 1; 00048 m_nbTimeBasisFunctions = -1; 00049 m2p_timeBasisFunctions = NULL; 00050 m2p_frameDurationsInMs = NULL; 00051 m2p_frameTimeStartInMs = NULL; 00052 m2p_frameTimeStopInMs = NULL; 00053 m_timeBasisFunctionsFile = ""; 00054 m_timeStaticFlag = false; 00055 m3p_quantificationFactors = NULL; 00056 m_ignoredCorrectionsList = ""; 00057 m_ignoreAttnCorrectionFlag = false; 00058 m_ignoreNormCorrectionFlag = false; 00059 m_ignoreRandCorrectionFlag = false; 00060 m_ignoreScatCorrectionFlag = false; 00061 m_ignoreDecaCorrectionFlag = false; 00062 m_ignoreBratCorrectionFlag = false; 00063 m_ignoreFdurCorrectionFlag = false; 00064 m_ignoreCaliCorrectionFlag = false; 00065 m_nbRespGates = 1; 00066 m_nbRespBasisFunctions = -1; 00067 m2p_respBasisFunctions = NULL; 00068 m_respBasisFunctionsFile = ""; 00069 m_respStaticFlag = false; 00070 m_nbCardGates = 1; 00071 m_nbCardBasisFunctions = -1; 00072 m2p_cardBasisFunctions = NULL; 00073 m_cardBasisFunctionsFile = ""; 00074 m_cardStaticFlag = false; 00075 m_nbThreadsForProjection = 1; 00076 m_nbThreadsForImageComputation = 1; 00077 m_mpiRank = 0; 00078 m_mpiSize = 1; 00079 m_nbBeds = -1; 00080 m_verbose = 0; 00081 m_checked = false; 00082 m_initialized = false; 00083 // Allocate Dynamic data manager object 00084 mp_DynamicDataManager = new oDynamicDataManager(); 00085 } 00086 00087 // ===================================================================== 00088 // --------------------------------------------------------------------- 00089 // --------------------------------------------------------------------- 00090 // ===================================================================== 00091 00092 oImageDimensionsAndQuantification::~oImageDimensionsAndQuantification() 00093 { 00094 // Free frame duration 00095 if (m2p_frameDurationsInMs) 00096 { 00097 for (int bed=0; bed<m_nbBeds; bed++) if (m2p_frameDurationsInMs[bed]) delete m2p_frameDurationsInMs[bed]; 00098 delete m2p_frameDurationsInMs; 00099 } 00100 // Free frame time start 00101 if (m2p_frameTimeStartInMs) 00102 { 00103 for (int bed=0; bed<m_nbBeds; bed++) if (m2p_frameTimeStartInMs[bed]) delete m2p_frameTimeStartInMs[bed]; 00104 delete m2p_frameTimeStartInMs; 00105 } 00106 // Free frame time stop 00107 if (m2p_frameTimeStopInMs) 00108 { 00109 for (int bed=0; bed<m_nbBeds; bed++) if (m2p_frameTimeStopInMs[bed]) delete m2p_frameTimeStopInMs[bed]; 00110 delete m2p_frameTimeStopInMs; 00111 } 00112 // Free quantification factors 00113 if (m3p_quantificationFactors) 00114 { 00115 for (int bed=0; bed<m_nbBeds; bed++) 00116 { 00117 if (m3p_quantificationFactors[bed]) 00118 { 00119 for (int fr=0; fr<m_nbTimeFrames; fr++) 00120 if (m3p_quantificationFactors[bed][fr]) delete m3p_quantificationFactors[bed][fr]; 00121 } 00122 delete m3p_quantificationFactors[bed]; 00123 } 00124 delete m3p_quantificationFactors; 00125 } 00126 // Delete the dynamic data manager 00127 delete mp_DynamicDataManager; 00128 } 00129 00130 // ===================================================================== 00131 // --------------------------------------------------------------------- 00132 // --------------------------------------------------------------------- 00133 // ===================================================================== 00134 00135 int oImageDimensionsAndQuantification::SetFlipOut(const string& a_flipOut) 00136 { 00137 // If empty, then put all flags to false and return 00138 if (a_flipOut=="") 00139 { 00140 m_flipOutX = false; 00141 m_flipOutY = false; 00142 m_flipOutZ = false; 00143 return 0; 00144 } 00145 // Test all possible settings ! 00146 if (a_flipOut=="x" || a_flipOut=="X") 00147 { 00148 m_flipOutX = true; 00149 m_flipOutY = false; 00150 m_flipOutZ = false; 00151 } 00152 else if (a_flipOut=="y" || a_flipOut=="Y") 00153 { 00154 m_flipOutX = false; 00155 m_flipOutY = true; 00156 m_flipOutZ = false; 00157 } 00158 else if (a_flipOut=="z" || a_flipOut=="Z") 00159 { 00160 m_flipOutX = false; 00161 m_flipOutY = false; 00162 m_flipOutZ = true; 00163 } 00164 else if (a_flipOut=="xy" || a_flipOut=="yx" || a_flipOut=="XY" || a_flipOut=="YX") 00165 { 00166 m_flipOutX = true; 00167 m_flipOutY = true; 00168 m_flipOutZ = false; 00169 } 00170 else if (a_flipOut=="zy" || a_flipOut=="yz" || a_flipOut=="ZY" || a_flipOut=="YZ") 00171 { 00172 m_flipOutX = false; 00173 m_flipOutY = true; 00174 m_flipOutZ = true; 00175 } 00176 else if (a_flipOut=="xz" || a_flipOut=="zx" || a_flipOut=="XZ" || a_flipOut=="ZX") 00177 { 00178 m_flipOutX = true; 00179 m_flipOutY = false; 00180 m_flipOutZ = true; 00181 } 00182 else if ( a_flipOut=="xyz" || a_flipOut=="xzy" || a_flipOut=="yxz" || a_flipOut=="yzx" || a_flipOut=="zxy" || a_flipOut=="zyx" || 00183 a_flipOut=="XYZ" || a_flipOut=="XZY" || a_flipOut=="YXZ" || a_flipOut=="YZX" || a_flipOut=="ZXY" || a_flipOut=="ZYX" ) 00184 { 00185 m_flipOutX = true; 00186 m_flipOutY = true; 00187 m_flipOutZ = true; 00188 } 00189 // Otherwise, throw an error 00190 else 00191 { 00192 Cerr("***** oImageDimensionsAndQuantification::SetFlipOut() -> Output flip settings is incorrect !" << endl); 00193 return 1; 00194 } 00195 // Normal end 00196 return 0; 00197 } 00198 00199 // ===================================================================== 00200 // --------------------------------------------------------------------- 00201 // --------------------------------------------------------------------- 00202 // ===================================================================== 00203 00204 int oImageDimensionsAndQuantification::SetNbThreads(const string& a_nbThreads) 00205 { 00206 // The number of threads can be given as two separated numbers to distinguish between the number of threads for projection computation 00207 // and the number of threads for image computation. Be careful that the thread-safe images are allocated with respect to the number 00208 // of projection threads; the number of threads used for image computation are strictly used for computation over images (i.e. when 00209 // using a loop over voxels. 00210 00211 // Look for a comma and check that there is only one comma without missing parameters 00212 size_t first_comma = a_nbThreads.find_first_of(","); 00213 size_t last_comma = a_nbThreads.find_last_of(","); 00214 if (first_comma!=last_comma || first_comma==0 || first_comma==a_nbThreads.length()-1) 00215 { 00216 Cerr("***** oImageDimensionsAndQuantification::SetNbThreads() -> Wrong syntax in the thread parameters ! See help." << endl); 00217 return 1; 00218 } 00219 00220 // Case for a single parameter 00221 if (first_comma==string::npos) 00222 { 00223 m_nbThreadsForProjection = atoi( a_nbThreads.c_str() ); 00224 m_nbThreadsForImageComputation = m_nbThreadsForProjection; 00225 } 00226 // Case for two parameters 00227 else 00228 { 00229 m_nbThreadsForProjection = atoi( (a_nbThreads.substr(0,first_comma)).c_str() ); 00230 m_nbThreadsForImageComputation = atoi( (a_nbThreads.substr(first_comma+1)).c_str() ); 00231 } 00232 00233 // Checks for negative threads numbers 00234 if (m_nbThreadsForProjection<0) 00235 { 00236 Cerr("***** oImageDimensionsAndQuantification::SetNbThreads() -> Negative number of threads provided for projection computation !" << endl); 00237 return 1; 00238 } 00239 if (m_nbThreadsForImageComputation<0) 00240 { 00241 Cerr("***** oImageDimensionsAndQuantification::SetNbThreads() -> Negative number of threads provided for image computation !" << endl); 00242 return 1; 00243 } 00244 00245 #ifdef CASTOR_OMP 00246 // If number of threads is 0, we automatically determine the maximum number of threads using OMP functions 00247 if (m_nbThreadsForProjection==0) m_nbThreadsForProjection = omp_get_max_threads(); 00248 if (m_nbThreadsForImageComputation==0) m_nbThreadsForImageComputation = omp_get_max_threads(); 00249 // At this step, set by default the number of threads for image operations 00250 omp_set_num_threads(m_nbThreadsForImageComputation); 00251 #endif 00252 00253 // End 00254 return 0; 00255 } 00256 00257 // ===================================================================== 00258 // --------------------------------------------------------------------- 00259 // --------------------------------------------------------------------- 00260 // ===================================================================== 00261 00262 int oImageDimensionsAndQuantification::CheckParameters() 00263 { 00264 // Number of threads 00265 if (m_nbThreadsForProjection<=0) 00266 { 00267 Cerr("***** oImageDimensionsAndQuantification::CheckParameters() -> Should provide a strictly positive number of threads !" << endl); 00268 return 1; 00269 } 00270 // Number of voxels 00271 if (m_nbVoxX<=0 || m_nbVoxY<=0 || m_nbVoxZ<=0) 00272 { 00273 Cerr("***** oImageDimensionsAndQuantification::CheckParameters() -> Should provide strictly positive number of voxels !" << endl); 00274 return 1; 00275 } 00276 // When FOV and voxel sizes are both not provided 00277 if ((m_voxSizeX<=0. || m_voxSizeY<=0. || m_voxSizeZ<=0.) && (m_fovSizeX<=0. || m_fovSizeY<=0. || m_fovSizeZ<=0.)) 00278 { 00279 Cerr("***** oImageDimensionsAndQuantification::CheckParameters() -> Should provide strictly positive voxel or FOV dimensions !" << endl); 00280 return 1; 00281 } 00282 // When both FOV and voxel sizes are provided 00283 if (m_voxSizeX>0. && m_voxSizeY>0. && m_voxSizeZ>0. && m_fovSizeX>0. && m_fovSizeY>0. && m_fovSizeZ>0.) 00284 { 00285 Cerr("***** oImageDimensionsAndQuantification::CheckParameters() -> Both FOV and voxels dimensions provided, should not provide both !" << endl); 00286 return 1; 00287 } 00288 // Check output transaxial FOV which is in percent (0 means we do not apply any FOV masking) 00289 if (m_fovOutPercent<0. || m_fovOutPercent>100.) 00290 { 00291 Cerr("***** oImageDimensionsAndQuantification::CheckParameters() -> Output transaxial FOV percentage must be in ]0:100] !" << endl); 00292 return 1; 00293 } 00294 // Check that the number of axial slices to be masked is between 0 and dimZ/2 00295 if (m_nbSliceOutMask<0 || m_nbSliceOutMask>m_nbVoxZ/2) 00296 { 00297 Cerr("***** oImageDimensionsAndQuantification::CheckParameters() -> Number of output axial slices to be masked is incorrectly set !" << endl); 00298 return 1; 00299 } 00300 // Check whether the DynamicDataManager object has been instanciated or not 00301 if (!mp_DynamicDataManager) 00302 { 00303 Cerr("***** oImageDimensionsAndQuantification::CheckParameters() -> Error : DynamicDataManager object not initialized !" << endl); 00304 return 1; 00305 } 00306 // Unlock 00307 m_checked = true; 00308 // End 00309 return 0; 00310 } 00311 00312 // ===================================================================== 00313 // --------------------------------------------------------------------- 00314 // --------------------------------------------------------------------- 00315 // ===================================================================== 00316 00317 int oImageDimensionsAndQuantification::Initialize() 00318 { 00319 // Mandatory check 00320 if (!m_checked) 00321 { 00322 Cerr("***** oImageDimensionsAndQuantification::Initialize() -> Cannot initialize before a call to CheckParameters() !" << endl); 00323 return 1; 00324 } 00325 00326 // Verbose 00327 if (m_verbose>=2) Cout("oImageDimensionsAndQuantification::Initialize() -> Initialize image dimensions, basis functions and quantification" << endl); 00328 00329 // Precompute number of voxels in a slice (XY) and total number of voxels (XYZ) 00330 m_nbVoxXY = m_nbVoxX * m_nbVoxY; 00331 m_nbVoxXYZ = m_nbVoxXY * m_nbVoxZ; 00332 00333 // If FOV dimensions are provided, then compute the voxel dimensions 00334 if (m_fovSizeX>0. && m_fovSizeY>0. && m_fovSizeZ>0.) 00335 { 00336 m_voxSizeX = m_fovSizeX / ((FLTNB)m_nbVoxX); 00337 m_voxSizeY = m_fovSizeY / ((FLTNB)m_nbVoxY); 00338 m_voxSizeZ = m_fovSizeZ / ((FLTNB)m_nbVoxZ); 00339 } 00340 // If voxel dimensions are provided, then compute the voxel dimensions 00341 else if (m_voxSizeX>0. && m_voxSizeY>0. && m_voxSizeZ>0.) 00342 { 00343 m_fovSizeX = m_voxSizeX * ((FLTNB)m_nbVoxX); 00344 m_fovSizeY = m_voxSizeY * ((FLTNB)m_nbVoxY); 00345 m_fovSizeZ = m_voxSizeZ * ((FLTNB)m_nbVoxZ); 00346 } 00347 00348 // Deal with frame list and quantification factors 00349 if (InitializeFramingAndQuantification()) 00350 { 00351 Cerr("***** oImageDimensionsAndQuantification::Initialize() -> A problem occured while initializing framing and quantification tabs !" << endl); 00352 return 1; 00353 } 00354 00355 // Deal with time basis functions 00356 if (InitializeTimeBasisFunctions()) 00357 { 00358 Cerr("***** oImageDimensionsAndQuantification::Initialize() -> A problem occured while initializing time basis functions !" << endl); 00359 return 1; 00360 } 00361 00362 // Deal with respiratory basis functions 00363 if (InitializeRespBasisFunctions()) 00364 { 00365 Cerr("***** oImageDimensionsAndQuantification::Initialize() -> A problem occured while initializing respiratory basis functions !" << endl); 00366 return 1; 00367 } 00368 00369 // Deal with cardiac basis functions 00370 if (InitializeCardBasisFunctions()) 00371 { 00372 Cerr("***** oImageDimensionsAndQuantification::Initialize() -> A problem occured while initializing cardiac basis functions !" << endl); 00373 return 1; 00374 } 00375 00376 // Deal with ignored corrections 00377 if (InitializeIgnoredCorrections()) 00378 { 00379 Cerr("***** oImageDimensionsAndQuantification::Initialize() -> A problem occured while initializing ignored corrections !" << endl); 00380 return 1; 00381 } 00382 00383 // Verbose 00384 if (m_verbose>=2) 00385 { 00386 Cout(" --> Image dimensions: [" << m_nbVoxX << ";" << m_nbVoxY << ";" << m_nbVoxZ << "] voxels of [" << m_voxSizeX << ";" << m_voxSizeY << ";" << m_voxSizeZ << "] mm3" << endl); 00387 Cout(" --> FOV size: [" << m_fovSizeX << ";" << m_fovSizeY << ";" << m_fovSizeZ << "] mm3" << endl); 00388 if (m_nbTimeFrames>1) 00389 { 00390 if (m_timeStaticFlag) Cout(" --> Time frames: " << m_nbTimeFrames << endl); 00391 else Cout(" --> Time frames: " << m_nbTimeFrames << " | Time basis functions: " << m_nbTimeBasisFunctions << endl); 00392 } 00393 if (m_respBasisFunctionsFile!="") 00394 { 00395 if (m_respStaticFlag) Cout(" --> Respiratory gates: " << m_nbRespGates << endl); 00396 else Cout(" --> Respiratory gates: " << m_nbRespGates << " | Respiratory basis functions: " << m_nbRespBasisFunctions << endl); 00397 } 00398 if (m_cardBasisFunctionsFile!="") 00399 { 00400 if (m_cardStaticFlag) Cout(" --> Cardiac gates: " << m_nbCardGates << endl); 00401 else Cout(" --> Cardiac gates: " << m_nbCardGates << " | Cardiac basis functions: " << m_nbCardBasisFunctions << endl); 00402 } 00403 if (m_nbThreadsForProjection>1 || m_nbThreadsForImageComputation>1) 00404 { 00405 if (m_nbThreadsForImageComputation==m_nbThreadsForProjection) Cout(" --> Number of parallel threads: " << m_nbThreadsForProjection << endl); 00406 else Cout(" --> Number of parallel threads for projection / image computation: [" << m_nbThreadsForProjection << "/" << m_nbThreadsForImageComputation << "]" << endl); 00407 } 00408 } 00409 00410 // Initialized 00411 m_initialized = true; 00412 00413 // End 00414 return 0; 00415 } 00416 00417 // ===================================================================== 00418 // --------------------------------------------------------------------- 00419 // --------------------------------------------------------------------- 00420 // ===================================================================== 00421 00422 int oImageDimensionsAndQuantification::InitializeFramingAndQuantification() 00423 { 00424 // TODO: for whole body 4D dynamic acquisitions, will need to read a set of frames for each bed where the number of frames and frame durations 00425 // should be the same for all beds, except which frame is dead and time start and stop also... Have to think about and manage this. 00426 // This function requires that the number of beds be set 00427 if (m_nbBeds<=0) 00428 { 00429 Cerr("***** oImageDimensionsAndQuantification::InitializeFramingAndQuantification() -> Number of beds must be set before setting the framing of the acquisition !" << endl); 00430 return 1; 00431 } 00432 00433 // -------------------------------------------------- 00434 // Allocate time start, stop and duration, to 1 frame 00435 m2p_frameDurationsInMs = (uint32_t**)malloc(m_nbBeds*sizeof(uint32_t*)); 00436 m2p_frameTimeStartInMs = (uint32_t**)malloc(m_nbBeds*sizeof(uint32_t*)); 00437 m2p_frameTimeStopInMs = (uint32_t**)malloc(m_nbBeds*sizeof(uint32_t*)); 00438 m3p_quantificationFactors = (FLTNB***)malloc(m_nbBeds*sizeof(FLTNB**)); 00439 for (int bed=0; bed<m_nbBeds; bed++) 00440 { 00441 m2p_frameDurationsInMs[bed] = (uint32_t*)malloc(1*sizeof(uint32_t)); 00442 m2p_frameTimeStartInMs[bed] = (uint32_t*)malloc(1*sizeof(uint32_t)); 00443 m2p_frameTimeStopInMs[bed] = (uint32_t*)malloc(1*sizeof(uint32_t)); 00444 } 00445 00446 // --------------------------------------- 00447 // Particular case where the list is empty 00448 if (m_frameList=="") 00449 { 00450 // Set number of frames to 1 00451 m_nbTimeFrames = 1; 00452 // Init them to 0 (the vDataFile will set them to the appropriate value later through the oImageDimensionsAndQuantification::SetAcquisitionTime() function) 00453 for (int bed=0; bed<m_nbBeds; bed++) 00454 { 00455 m2p_frameDurationsInMs[bed][0] = 0; 00456 m2p_frameTimeStartInMs[bed][0] = 0; 00457 m2p_frameTimeStopInMs[bed][0] = 0; 00458 // Allocate quantification factors and set them to 1. 00459 m3p_quantificationFactors[bed] = (FLTNB**)malloc(1*sizeof(FLTNB*)); 00460 m3p_quantificationFactors[bed][0] = (FLTNB*)malloc(m_nbRespGates*m_nbCardGates*sizeof(FLTNB)); 00461 for (int g=0; g<m_nbRespGates*m_nbCardGates; g++) m3p_quantificationFactors[bed][0][g] = 1.; 00462 } 00463 // Exit the function 00464 return 0; 00465 } 00466 00467 // -------------------------------------------------------------------- 00468 // Otherwise, get the parameter as a non-constant string and process it 00469 string frame_list = m_frameList; 00470 // Current time start 00471 uint32_t time_start = 0; 00472 // Init number of frames 00473 m_nbTimeFrames = 0; 00474 // Loop on all commas found in the list 00475 size_t comma_pos = 0; 00476 while ((comma_pos=frame_list.find_first_of(","))!=string::npos) 00477 { 00478 // Extract the parameter before the comma 00479 string param = frame_list.substr(0,comma_pos); 00480 // Search for a dead-frame sign ";0" 00481 size_t semi_column_pos = param.find(";0"); 00482 // If does not find then this is an actual frame 00483 if (semi_column_pos==string::npos) 00484 { 00485 // Increment number of frames 00486 m_nbTimeFrames++; 00487 // Realloc the time start, stop and duration 00488 for (int bed=0; bed<m_nbBeds; bed++) 00489 { 00490 m2p_frameDurationsInMs[bed] = (uint32_t*)realloc( m2p_frameDurationsInMs[bed], m_nbTimeFrames*sizeof(uint32_t) ); 00491 m2p_frameTimeStartInMs[bed] = (uint32_t*)realloc( m2p_frameTimeStartInMs[bed], m_nbTimeFrames*sizeof(uint32_t) ); 00492 m2p_frameTimeStopInMs[bed] = (uint32_t*)realloc( m2p_frameTimeStopInMs[bed] , m_nbTimeFrames*sizeof(uint32_t) ); 00493 } 00494 // Extract frame duration 00495 FLTNB duration = atof(param.c_str()); 00496 if (duration<=0.) 00497 { 00498 Cerr("***** oImageDimensionsAndQuantification::InitializeFramingAndQuantification() -> Found a duration which is negative or null !" << endl); 00499 return 1; 00500 } 00501 // Affect time start, stop and duration 00502 for (int bed=0; bed<m_nbBeds; bed++) 00503 { 00504 m2p_frameDurationsInMs[bed][m_nbTimeFrames-1] = ((uint32_t)(duration*1000.)); 00505 m2p_frameTimeStartInMs[bed][m_nbTimeFrames-1] = time_start; 00506 } 00507 time_start += m2p_frameDurationsInMs[0][m_nbTimeFrames-1]; 00508 for (int bed=0; bed<m_nbBeds; bed++) 00509 m2p_frameTimeStopInMs[bed][m_nbTimeFrames-1] = time_start; 00510 } 00511 // Otherwise, this is a dead frame 00512 else 00513 { 00514 // Check that the keyword ";0" if at the end of the parameter 00515 if (semi_column_pos!=param.size()-2) 00516 { 00517 Cerr("***** oImageDimensionsAndQuantification::InitializeFramingAndQuantification() -> Invalid dead frame format '" << param << "' !" << endl); 00518 return 1; 00519 } 00520 // Extract the dead frame duration 00521 FLTNB duration = atof( (param.substr(0,semi_column_pos)).c_str() ); 00522 if (duration<=0.) 00523 { 00524 Cerr("***** oImageDimensionsAndQuantification::InitializeFramingAndQuantification() -> Found a duration which is negative or null !" << endl); 00525 return 1; 00526 } 00527 // Increment time start 00528 time_start += ((uint32_t)(duration*1000.)); 00529 } 00530 // Remove first parameter 00531 frame_list = frame_list.substr(comma_pos+1); 00532 } 00533 // Last parameter to extract; search for a dead-frame sign ";0" 00534 size_t semi_column_pos = frame_list.find(";0"); 00535 // If does not find then this is an actual frame 00536 if (semi_column_pos==string::npos) 00537 { 00538 // Increment number of frames 00539 m_nbTimeFrames++; 00540 // Realloc the time start, stop and duration 00541 for (int bed=0; bed<m_nbBeds; bed++) 00542 { 00543 m2p_frameDurationsInMs[bed] = (uint32_t*)realloc(m2p_frameDurationsInMs[bed],m_nbTimeFrames*sizeof(uint32_t)); 00544 m2p_frameTimeStartInMs[bed] = (uint32_t*)realloc(m2p_frameTimeStartInMs[bed],m_nbTimeFrames*sizeof(uint32_t)); 00545 m2p_frameTimeStopInMs[bed] = (uint32_t*)realloc(m2p_frameTimeStopInMs[bed] ,m_nbTimeFrames*sizeof(uint32_t)); 00546 } 00547 // Extract frame duration 00548 FLTNB duration = atof(frame_list.c_str()); 00549 if (duration<=0.) 00550 { 00551 Cerr("***** oImageDimensionsAndQuantification::InitializeFramingAndQuantification() -> Found a duration which is negative or null !" << endl); 00552 return 1; 00553 } 00554 // Affect time start, stop and duration 00555 for (int bed=0; bed<m_nbBeds; bed++) 00556 { 00557 m2p_frameDurationsInMs[bed][m_nbTimeFrames-1] = ((uint32_t)(duration*1000.)); 00558 m2p_frameTimeStartInMs[bed][m_nbTimeFrames-1] = time_start; 00559 } 00560 time_start += m2p_frameDurationsInMs[0][m_nbTimeFrames-1]; 00561 for (int bed=0; bed<m_nbBeds; bed++) 00562 m2p_frameTimeStopInMs[bed][m_nbTimeFrames-1] = time_start; 00563 } 00564 // Otherwise, this is a dead frame, so we just ignore it as it is in last position 00565 else 00566 { 00567 ; 00568 } 00569 00570 // ---------------------------------------------- 00571 // Allocate and affect the quantification factors 00572 for (int bed=0; bed<m_nbBeds; bed++) 00573 { 00574 m3p_quantificationFactors[bed] = (FLTNB**)malloc(m_nbTimeFrames*sizeof(FLTNB*)); 00575 for (int fr=0; fr<m_nbTimeFrames; fr++) 00576 { 00577 m3p_quantificationFactors[bed][fr] = (FLTNB*)malloc(m_nbRespGates*m_nbCardGates*sizeof(FLTNB)); 00578 for (int g=0; g<m_nbRespGates*m_nbCardGates; g++) m3p_quantificationFactors[bed][fr][g] = 1.; 00579 } 00580 } 00581 00582 // --- 00583 // End 00584 return 0; 00585 } 00586 00587 // ===================================================================== 00588 // --------------------------------------------------------------------- 00589 // --------------------------------------------------------------------- 00590 // ===================================================================== 00591 00592 int oImageDimensionsAndQuantification::InitializeTimeBasisFunctions() 00593 { 00594 // Case 1: no time basis functions file provided 00595 if (m_timeBasisFunctionsFile=="") 00596 { 00597 // Set the number of time basis functions to the number of frames 00598 m_nbTimeBasisFunctions = m_nbTimeFrames; 00599 // Allocate the conversion matrix from basis functions to frames 00600 m2p_timeBasisFunctions = (FLTNB**)malloc(m_nbTimeBasisFunctions*sizeof(FLTNB*)); 00601 for (int tbf=0; tbf<m_nbTimeBasisFunctions; tbf++) 00602 { 00603 m2p_timeBasisFunctions[tbf] = (FLTNB*)malloc(m_nbTimeFrames*sizeof(FLTNB)); 00604 for (int fr=0; fr<m_nbTimeFrames; fr++) 00605 { 00606 // Set diagonal to 1, the rest to 0 00607 if (tbf==fr) m2p_timeBasisFunctions[tbf][fr] = 1.; 00608 else m2p_timeBasisFunctions[tbf][fr] = 0.; 00609 } 00610 } 00611 // Set the static flag to true 00612 m_timeStaticFlag = true; 00613 } 00614 // Case 2: there is a file 00615 else 00616 { 00617 // Allocate the conversion matrix from basis functions to frames 00618 m2p_timeBasisFunctions = (FLTNB**)malloc(m_nbTimeBasisFunctions*sizeof(FLTNB*)); 00619 for (int tbf=0; tbf<m_nbTimeBasisFunctions; tbf++) m2p_timeBasisFunctions[tbf] = (FLTNB*)malloc(m_nbTimeFrames*sizeof(FLTNB)); 00620 // Open file 00621 ifstream fin(m_timeBasisFunctionsFile.c_str()); 00622 if (!fin) 00623 { 00624 Cerr("***** oImageDimensionsAndQuantification::InitializeTimeBasisFunctions() -> Input time basis functions file '" << m_timeBasisFunctionsFile << "' is missing or corrupted !" << endl); 00625 return 1; 00626 } 00627 // Read values 00628 for (int tbf=0; tbf<m_nbTimeBasisFunctions; tbf++) for (int fr=0; fr<m_nbTimeFrames; fr++) fin >> m2p_timeBasisFunctions[tbf][fr]; 00629 // Close file 00630 fin.close(); 00631 } 00632 // End 00633 return 0; 00634 } 00635 00636 // ===================================================================== 00637 // --------------------------------------------------------------------- 00638 // --------------------------------------------------------------------- 00639 // ===================================================================== 00640 00641 int oImageDimensionsAndQuantification::InitializeRespBasisFunctions() 00642 { 00643 // Case 1: no respiratory basis functions file provided 00644 if (m_respBasisFunctionsFile=="") 00645 { 00646 // Set the number of respiratory basis functions to the number of respiratory gates 00647 m_nbRespBasisFunctions = m_nbRespGates; 00648 00649 // Allocate the conversion matrix from basis functions to respiratory gates 00650 m2p_respBasisFunctions = (FLTNB**)malloc(m_nbRespBasisFunctions*sizeof(FLTNB*)); 00651 for (int rbf=0; rbf<m_nbRespBasisFunctions; rbf++) 00652 { 00653 m2p_respBasisFunctions[rbf] = (FLTNB*)malloc(m_nbRespGates*sizeof(FLTNB)); 00654 for (int rg=0; rg<m_nbRespGates; rg++) 00655 { 00656 // Set diagonal to 1, the rest to 0 00657 if (rbf==rg) m2p_respBasisFunctions[rbf][rg] = 1.; 00658 else m2p_respBasisFunctions[rbf][rg] = 0.; 00659 } 00660 } 00661 // Set the static flag to true 00662 m_respStaticFlag = true; 00663 } 00664 // Case 2: there is a file 00665 else 00666 { 00667 // Allocate the conversion matrix from basis functions to respiratory gates 00668 m2p_respBasisFunctions = (FLTNB**)malloc(m_nbRespBasisFunctions*sizeof(FLTNB*)); 00669 for (int rbf=0; rbf<m_nbRespBasisFunctions; rbf++) m2p_respBasisFunctions[rbf] = (FLTNB*)malloc(m_nbRespGates*sizeof(FLTNB)); 00670 // Open file 00671 ifstream fin(m_respBasisFunctionsFile.c_str()); 00672 if (!fin) 00673 { 00674 Cerr("***** oImageDimensionsAndQuantification::InitializeRespBasisFunctions() -> Input respiratory basis functions file '" << m_respBasisFunctionsFile << "' is missing or corrupted !" << endl); 00675 return 1; 00676 } 00677 // Read values 00678 for (int rbf=0; rbf<m_nbRespBasisFunctions; rbf++) for (int rg=0; rg<m_nbRespGates; rg++) fin >> m2p_respBasisFunctions[rbf][rg]; 00679 // Close file 00680 fin.close(); 00681 } 00682 // End 00683 return 0; 00684 } 00685 00686 // ===================================================================== 00687 // --------------------------------------------------------------------- 00688 // --------------------------------------------------------------------- 00689 // ===================================================================== 00690 00691 int oImageDimensionsAndQuantification::InitializeCardBasisFunctions() 00692 { 00693 // Case 1: no cardiac basis functions file provided 00694 if (m_cardBasisFunctionsFile=="") 00695 { 00696 // Set the number of cardiac basis functions to the number of cardiac gates 00697 m_nbCardBasisFunctions = m_nbCardGates; 00698 // Allocate the conversion matrix from cardiac basis functions to cardiac gates 00699 m2p_cardBasisFunctions = (FLTNB**)malloc(m_nbCardBasisFunctions*sizeof(FLTNB*)); 00700 for (int cbf=0; cbf<m_nbCardBasisFunctions; cbf++) 00701 { 00702 m2p_cardBasisFunctions[cbf] = (FLTNB*)malloc(m_nbCardGates*sizeof(FLTNB)); 00703 for (int cg=0; cg<m_nbCardGates; cg++) 00704 { 00705 // Set diagonal to 1, the rest to 0 00706 if (cbf==cg) m2p_cardBasisFunctions[cbf][cg] = 1.; 00707 else m2p_cardBasisFunctions[cbf][cg] = 0.; 00708 } 00709 } 00710 // Set the static flag to true 00711 m_cardStaticFlag = true; 00712 } 00713 // Case 2: there is a file 00714 else 00715 { 00716 // Allocate the conversion matrix from basis functions to cardiac gates 00717 m2p_cardBasisFunctions = (FLTNB**)malloc(m_nbCardBasisFunctions*sizeof(FLTNB*)); 00718 for (int cbf=0; cbf<m_nbCardBasisFunctions; cbf++) m2p_cardBasisFunctions[cbf] = (FLTNB*)malloc(m_nbCardGates*sizeof(FLTNB)); 00719 // Open file 00720 ifstream fin(m_cardBasisFunctionsFile.c_str()); 00721 if (!fin) 00722 { 00723 Cerr("***** oImageDimensionsAndQuantification::InitializeCardBasisFunctions() -> Input cardiac basis functions file '" << m_cardBasisFunctionsFile << "' is missing or corrupted !" << endl); 00724 return 1; 00725 } 00726 // Read values 00727 for (int cbf=0; cbf<m_nbCardBasisFunctions; cbf++) for (int cg=0; cg<m_nbCardGates; cg++) fin >> m2p_cardBasisFunctions[cbf][cg]; 00728 // Close file 00729 fin.close(); 00730 } 00731 // End 00732 return 0; 00733 } 00734 00735 // ===================================================================== 00736 // --------------------------------------------------------------------- 00737 // --------------------------------------------------------------------- 00738 // ===================================================================== 00739 00740 int oImageDimensionsAndQuantification::InitializeIgnoredCorrections() 00741 { 00742 // If the m_ignoredCorrectionsList is empty, force all booleans to false and leave this function 00743 if (m_ignoredCorrectionsList=="") 00744 { 00745 m_ignoreAttnCorrectionFlag = false; 00746 m_ignoreNormCorrectionFlag = false; 00747 m_ignoreRandCorrectionFlag = false; 00748 m_ignoreScatCorrectionFlag = false; 00749 m_ignoreDecaCorrectionFlag = false; 00750 m_ignoreBratCorrectionFlag = false; 00751 m_ignoreFdurCorrectionFlag = false; 00752 m_ignoreCaliCorrectionFlag = false; 00753 return 0; 00754 } 00755 00756 // Then, count the number of commas in the m_ignoredCorrectionsList to know the number of keywords 00757 size_t nb_keywords = count(m_ignoredCorrectionsList.begin(), m_ignoredCorrectionsList.end(), ',') + 1; 00758 00759 // Read the keywords 00760 string *p_keywords = new string[nb_keywords]; 00761 00762 if (ReadStringOption(m_ignoredCorrectionsList, p_keywords, nb_keywords, ",", "Ignored corrections list")) 00763 { 00764 Cerr("***** oImageDimensionsAndQuantification::InitializeIgnoredCorrections() -> An error occured while reading the list of ignored corrections !" << endl); 00765 return 1; 00766 } 00767 00768 // Process them 00769 for (size_t k=0; k<nb_keywords; k++) 00770 { 00771 // Test the different keywords 00772 if (p_keywords[k]=="attn") m_ignoreAttnCorrectionFlag = true; 00773 else if (p_keywords[k]=="norm") m_ignoreNormCorrectionFlag = true; 00774 else if (p_keywords[k]=="scat") m_ignoreScatCorrectionFlag = true; 00775 else if (p_keywords[k]=="rand") m_ignoreRandCorrectionFlag = true; 00776 else if (p_keywords[k]=="deca") m_ignoreDecaCorrectionFlag = true; 00777 else if (p_keywords[k]=="brat") m_ignoreBratCorrectionFlag = true; 00778 else if (p_keywords[k]=="fdur") m_ignoreFdurCorrectionFlag = true; 00779 else if (p_keywords[k]=="cali") m_ignoreCaliCorrectionFlag = true; 00780 else 00781 { 00782 Cerr("***** oImageDimensionsAndQuantification::InitializeIgnoredCorrections() -> Unknown keyword '" << p_keywords[k] << "' in the provided ignored corrections list !" << endl); 00783 return 1; 00784 } 00785 } 00786 00787 delete[] p_keywords; 00788 // Normal end 00789 return 0; 00790 } 00791 00792 // ===================================================================== 00793 // --------------------------------------------------------------------- 00794 // --------------------------------------------------------------------- 00795 // ===================================================================== 00796 00797 int oImageDimensionsAndQuantification::SetAcquisitionTime(int a_bed, FLTNB a_timeStartInSec, FLTNB a_durationInSec) 00798 { 00799 // Check if initialized 00800 if (!m_initialized) 00801 { 00802 Cerr("***** oImageDimensionsAndQuantification::SetAcquisitionTime() -> Object not initialized !" << endl); 00803 return 1; 00804 } 00805 // If the number of frames is more than one, then it was already specified, so we just update the quantification factors and exit 00806 if (m_nbTimeFrames>1) 00807 { 00808 if (!m_ignoreFdurCorrectionFlag) 00809 { 00810 for (int fr=0; fr<m_nbTimeFrames; fr++) for(int g=0 ; g<m_nbRespGates*m_nbCardGates ; g++) 00811 m3p_quantificationFactors[a_bed][fr][g] /= (((FLTNB)(m2p_frameDurationsInMs[a_bed][fr]))/1000.); 00812 } 00813 } 00814 // If timeStart, timeStop and duration are 0, then it means that it was not specified yet 00815 else if (m2p_frameDurationsInMs[a_bed][0]==0 && m2p_frameTimeStartInMs[a_bed][0]==0 && m2p_frameTimeStopInMs[a_bed][0]==0) 00816 { 00817 m2p_frameDurationsInMs[a_bed][0] = ((uint32_t)(a_durationInSec*1000.)); 00818 m2p_frameTimeStartInMs[a_bed][0] = ((uint32_t)(a_timeStartInSec*1000.)); 00819 m2p_frameTimeStopInMs[a_bed][0] = m2p_frameTimeStartInMs[a_bed][0] + m2p_frameDurationsInMs[a_bed][0]; 00820 // Apply frame's duration onto quantification factors 00821 if (!m_ignoreFdurCorrectionFlag) 00822 { 00823 for(int g=0 ; g<m_nbRespGates*m_nbCardGates ; g++) 00824 m3p_quantificationFactors[a_bed][0][g] /= a_durationInSec; 00825 } 00826 } 00827 // Verbose 00828 if (m_verbose>=2 && (a_bed==m_nbBeds-1)) 00829 { 00830 // Case 1: a static single bed 00831 if (m_nbTimeFrames==1 && m_nbBeds==1) 00832 { 00833 Cout("oImageDimensionsAndQuantification::SetAcquisitionTime() -> Static single bed acquisition with duration [ " << GetFrameTimeStartInSec(0,0) << " : " 00834 << GetFrameTimeStopInSec(0,0) << " ] seconds" << endl); 00835 } 00836 // Case 2: a static multi beds 00837 else if (m_nbTimeFrames==1 && m_nbBeds>1) 00838 { 00839 Cout("oImageDimensionsAndQuantification::SetAcquisitionTime() -> Static " << m_nbBeds << " beds acquisition with following bed durations:" << endl); 00840 for (int bed=0; bed<m_nbBeds; bed++) 00841 Cout(" --> Bed " << bed+1 << " with duration [ " << GetFrameTimeStartInSec(bed,0) << " : " << GetFrameTimeStopInSec(bed,0) << " ] seconds" << endl); 00842 } 00843 // Case 3: a dynamic single bed 00844 else if (m_nbTimeFrames>1 && m_nbBeds==1) 00845 { 00846 Cout("oImageDimensionsAndQuantification::SetAcquisitionTime() -> Dynamic single bed acquisition with following " << m_nbTimeFrames << " frame durations:" << endl); 00847 for (int fr=0; fr<m_nbTimeFrames; fr++) 00848 Cout(" --> Frame " << fr+1 << " with duration [ " << GetFrameTimeStartInSec(0,fr) << " : " << GetFrameTimeStopInSec(0,fr) << " ] seconds" << endl); 00849 } 00850 // Case 4: a dynamic multi beds 00851 else 00852 { 00853 Cout("oImageDimensionsAndQuantification::SetAcquisitionTime() -> Dynamic " << m_nbBeds << " beds acquistion with following " << m_nbTimeFrames << " frame durations:" << endl); 00854 for (int bed=0; bed<m_nbBeds; bed++) 00855 { 00856 Cout(" --> Bed " << bed+1 << " as following framing:" << endl); 00857 for (int fr=0; fr<m_nbTimeFrames; fr++) 00858 Cout(" | Frame " << fr+1 << " with duration [ " << GetFrameTimeStartInSec(bed,fr) << " : " << GetFrameTimeStopInSec(bed,fr) << " ] seconds" << endl); 00859 } 00860 } 00861 // Correct for frame duration or not 00862 if (m_ignoreFdurCorrectionFlag) Cout(" --> Ignore frame duration correction" << endl); 00863 else Cout(" --> Correct for frame duration" << endl); 00864 } 00865 // End 00866 return 0; 00867 } 00868 00869 // ===================================================================== 00870 // --------------------------------------------------------------------- 00871 // --------------------------------------------------------------------- 00872 // ===================================================================== 00873 00874 int oImageDimensionsAndQuantification::SetCalibrationFactor(int a_bed, FLTNB a_calibrationFactor) 00875 { 00876 // Check if initialized 00877 if (!m_initialized) 00878 { 00879 Cerr("***** oImageDimensionsAndQuantification::SetCalibrationFactor() -> Object not initialized !" << endl); 00880 return 1; 00881 } 00882 // Ignore calibration factor if asked for 00883 if (m_ignoreCaliCorrectionFlag) 00884 { 00885 // Verbose 00886 if (m_verbose>=2 && a_bed==m_nbBeds-1) Cout("oImageDimensionsAndQuantification::SetCalibrationFactor() -> Ignore calibration factor correction" << endl); 00887 // Exit the function 00888 return 0; 00889 } 00890 // Check if calibration factor is strictly positive 00891 if (a_calibrationFactor<=0.) 00892 { 00893 Cerr("***** oImageDimensionsAndQuantification::SetCalibrationFactor() -> Provided calibration factor (" << a_calibrationFactor << ") is negative or null !" << endl); 00894 return 1; 00895 } 00896 // Affect quantification factor for all frames for this bed (even though the calibration factor should be the same 00897 // for all beds, we do not check it, it should be self consistent in the input files) 00898 if (!m_ignoreCaliCorrectionFlag) 00899 { 00900 for (int fr=0; fr<m_nbTimeFrames; fr++) for(int g=0 ; g<m_nbRespGates*m_nbCardGates ; g++) 00901 m3p_quantificationFactors[a_bed][fr][g] *= a_calibrationFactor; 00902 } 00903 // Verbose 00904 if (m_verbose>=2 && a_bed==m_nbBeds-1) 00905 Cout("oImageDimensionsAndQuantification::SetCalibrationFactor() -> Correct for following calibration factor: " << a_calibrationFactor << endl); 00906 // End 00907 return 0; 00908 } 00909 00910 // ===================================================================== 00911 // --------------------------------------------------------------------- 00912 // --------------------------------------------------------------------- 00913 // ===================================================================== 00914 00915 int oImageDimensionsAndQuantification::SetDynamicSpecificQuantificationFactors(const string& a_quantificationFile) 00916 { 00917 // Check if initialized 00918 if (!m_initialized) 00919 { 00920 Cerr("***** oImageDimensionsAndQuantification::SetDynamicSpecificQuantificationFactors() -> Object not initialized !" << endl); 00921 return 1; 00922 } 00923 00924 // Provide the DynamicDataManager with a 2D array containing quantitative factors to update during initialization 00925 mp_DynamicDataManager->SetDynamicSpecificQuantificationFactors(m3p_quantificationFactors[0]); 00926 00927 // Exit the function if no file provided 00928 if (a_quantificationFile=="") return 0; 00929 00930 // Currently, it is not possible to reconstruction a multi-bed gated acquisition, so we throw an error in this function 00931 if (m_nbBeds>1) 00932 { 00933 Cerr("***** oImageDimensionsAndQuantification::SetDynamicSpecificQuantificationFactors() -> Multi-bed gated acquisitions cannot be reconstructed !" << endl); 00934 return 1; 00935 } 00936 00937 // Verbose 00938 if (m_verbose>=2) Cout("oImageDimensionsAndQuantification::SetDynamicSpecificQuantificationFactors()-> Processing quantification file '" << a_quantificationFile << "'" << endl); 00939 00940 // GET QUANTIFICATION CORRECTION FACTORS WITH RESPECT TO FRAMES or GATES 00941 00942 // READ USER-DEFINED QUANTIFICATION FACTORS 00943 00944 // TODO : For now we assume the configuration file holds : 00945 // TODO : - frame quantification factors : specific to each bed & frame 00946 // TODO : - gate quantification factors : specific to each frame and gates 00947 // TODO : Might require a new ReadDataASCIIFile function to recover this 00948 // 00949 // TODO : Additionally, we should be careful in documentation and explained which gate quantitfication factors are natively taken into account or not 00950 00951 00952 // Quantification factors for each bed/frame/gate, intialized to 1. 00953 FLTNB*** pp_dynamic_quantification_factors = new FLTNB**[m_nbBeds]; 00954 for (int bed=0 ; bed<m_nbBeds ; bed++) 00955 { 00956 pp_dynamic_quantification_factors[bed] = new FLTNB*[m_nbTimeFrames]; 00957 for (int fr=0; fr<m_nbTimeFrames; fr++) 00958 { 00959 pp_dynamic_quantification_factors[bed][fr] = new FLTNB[m_nbRespGates*m_nbCardGates]; 00960 for (int g=0; g<m_nbRespGates*m_nbCardGates; g++) pp_dynamic_quantification_factors[bed][fr][g] = 1.; 00961 } 00962 } 00963 00964 // Build bed-related key words 00965 string *bed_name = new string[m_nbBeds + 1]; 00966 for(int bed=0 ; bed<m_nbBeds ; bed++) 00967 { 00968 ostringstream oss( ostringstream::out ); 00969 oss << "bed" << bed+1; 00970 bed_name[bed] = oss.str(); 00971 } 00972 bed_name[m_nbBeds] = "eof"; 00973 00974 for (int bed=0 ; bed<m_nbBeds ; bed++) 00975 { 00976 int return_value = ReadDataASCIIFile(a_quantificationFile, "QUANTIFICATION_FACTORS", pp_dynamic_quantification_factors[bed], m_nbRespGates*m_nbCardGates, m_nbTimeFrames, KEYWORD_MANDATORY, bed_name[bed], bed_name[bed+1]); 00977 if (return_value<0) // string not found error 00978 { 00979 Cerr("***** oImageDimensionsAndQuantification::SetDynamicSpecificQuantificationFactors() -> Didn't found quantitative factors in file " << a_quantificationFile << " !" << endl); 00980 return 1; 00981 } 00982 else if(return_value == 1) // reading error 00983 { 00984 Cerr("***** oImageDimensionsAndQuantification::SetDynamicSpecificQuantificationFactors() -> An error occured while trying to recover specific quantitative factors for frame !" << endl); 00985 return 1; 00986 } 00987 else if(return_value == 0) // correct reading 00988 { 00989 for(int bed=0 ; bed<m_nbBeds ; bed++) 00990 for (int fr=0; fr<m_nbTimeFrames; fr++) 00991 for (int g=0; g<m_nbRespGates*m_nbCardGates; g++) 00992 if (pp_dynamic_quantification_factors[bed][fr][g] <= 0) 00993 { 00994 Cerr("***** oImageDimensionsAndQuantification::SetDynamicSpecificQuantificationFactors() -> Provided quantification factor (" << pp_dynamic_quantification_factors[bed][fr][g] << ") is negative or null !" << endl); 00995 return 1; 00996 } 00997 else 00998 { 00999 m3p_quantificationFactors[bed][fr][g] *= pp_dynamic_quantification_factors[bed][fr][g]; 01000 } 01001 } 01002 } 01003 01004 // Delete temporary tabs 01005 for (int bed=0; bed<m_nbBeds; bed++) 01006 { 01007 for (int fr=0; fr<m_nbTimeFrames; fr++) delete pp_dynamic_quantification_factors[bed][fr]; 01008 delete[] pp_dynamic_quantification_factors[bed]; 01009 } 01010 delete[] pp_dynamic_quantification_factors; 01011 delete bed_name; 01012 01013 // End 01014 return 0; 01015 } 01016 01017 // ===================================================================== 01018 // --------------------------------------------------------------------- 01019 // --------------------------------------------------------------------- 01020 // ===================================================================== 01021 01022 FLTNB oImageDimensionsAndQuantification::GetQuantificationFactor(int a_bed, int a_frame, int a_respGate, int a_cardGate) 01023 { 01024 // Make sure that if respiratory or cardiac motion is enabled, we will recover the factor corresponding to the first image 01025 // TODO TM : perhaps more logical to initialize the oImageDimensions with the actual number of gates even when motion is enabled 01026 // TODO TM : (i.e : m_nbRespGates or m_nbCardGates will be equal to 1, but the sensitivity image (list-mode) will require quantification factor for each gate 01027 // TODO TM : and add GetNbRespImages() functions anywhere in the code we need to know the number of resp/card images to be reconstructed 01028 // SS: If the motion correction is enabled, then all gates are pulled together in the optimization process, so there is no need for any quantification factor specific to each gate... 01029 if (m_nbRespGates == 1) a_respGate = 0; 01030 if (m_nbCardGates == 1) a_cardGate = 0; 01031 return m3p_quantificationFactors[a_bed][a_frame][a_respGate*m_nbCardGates+a_cardGate]; 01032 } 01033 01034 // ===================================================================== 01035 // --------------------------------------------------------------------- 01036 // --------------------------------------------------------------------- 01037 // ===================================================================== 01038 01039 int oImageDimensionsAndQuantification::SetSPECTIsotope(int a_bed, const string& a_isotope) 01040 { 01041 // Check if initialized 01042 if (!m_initialized) 01043 { 01044 Cerr("***** oImageDimensionsAndQuantification::SetSPECTIsotope() -> Object not initialized !" << endl); 01045 return 1; 01046 } 01047 01048 // Not yet implemented 01049 01050 // Normal end 01051 return 0; 01052 } 01053 01054 // ===================================================================== 01055 // --------------------------------------------------------------------- 01056 // --------------------------------------------------------------------- 01057 // ===================================================================== 01058 01059 int oImageDimensionsAndQuantification::SetPETIsotope(int a_bed, const string& a_isotope) 01060 { 01061 // Check if initialized 01062 if (!m_initialized) 01063 { 01064 Cerr("***** oImageDimensionsAndQuantification::SetPETIsotope() -> Object not initialized !" << endl); 01065 return 1; 01066 } 01067 01068 // ------------------------------------------------------------------ 01069 // Preliminary step is checking if we ignore isotope corrections 01070 // ------------------------------------------------------------------ 01071 01072 if (m_ignoreBratCorrectionFlag && m_ignoreDecaCorrectionFlag) 01073 { 01074 // Verbose 01075 if (m_verbose>=2 && a_bed==m_nbBeds-1) 01076 Cout("oImageDimensionsAndQuantification::SetPETIsotope() -> Ignore isotope dependent corrections" << endl); 01077 // Exit the function 01078 return 0; 01079 } 01080 01081 // ------------------------------------------------------------------ 01082 // Preliminary step is checking for unknown keyword 01083 // ------------------------------------------------------------------ 01084 01085 // Check if the isotope is named as unknown 01086 if (a_isotope=="UNKNOWN" || a_isotope=="Unknown" || a_isotope=="unknown") 01087 { 01088 // In this case we simply assume perfect branching ratio and infinite half-life 01089 if (m_verbose>=2 && a_bed==m_nbBeds-1) 01090 Cout("oImageDimensionsAndQuantification::SetPETIsotope() -> Un-specified isotope; no decay nor branching ratio correction" << endl); 01091 // Exit the function 01092 return 0; 01093 } 01094 01095 // ------------------------------------------------------------------ 01096 // First step is open the isotopes file and find the provided isotope 01097 // ------------------------------------------------------------------ 01098 01099 // Get oOutputManager instance and config directory 01100 sOutputManager* p_output_manager = sOutputManager::GetInstance(); 01101 string config_dir = p_output_manager->GetPathToConfigDir(); 01102 01103 // Open the isotope file based on config directory 01104 string file_name = config_dir + "/misc/isotopes_pet.txt"; 01105 ifstream fin(file_name.c_str()); 01106 01107 if (!fin) 01108 { 01109 Cerr("***** oImageDimensionsAndQuantification::SetPETIsotope() -> Failed to open PET isotopes data file '" << file_name << "' !" << endl); 01110 return 1; 01111 } 01112 01113 // Loop on lines to find the isotope 01114 int line_max_size = 10240; 01115 char *line = new char[line_max_size]; 01116 FLTNB half_life = -1.; 01117 FLTNB branching_ratio = -1.; 01118 bool found_it = false; 01119 fin.getline(line,line_max_size); 01120 while (!fin.eof()) 01121 { 01122 // For the word position 01123 size_t found_position; 01124 // Cast the line into a string for convenience 01125 string test = (string)line; 01126 // Jump to next line if we find a # character as the first character 01127 if ((found_position=test.find("#"))==0) 01128 { 01129 // Read a new line before continuing 01130 fin.getline(line,line_max_size); 01131 continue; 01132 } 01133 // Check if we see the isotope name 01134 found_position = test.find(a_isotope); 01135 if (found_position!=string::npos) 01136 { 01137 // Each line is organised as follows: (some spaces/tabs) isotope_name (some spaces/tabs) half_life (some spaces/tabs) branching_ratio 01138 // We first remove the isotope name from the line 01139 test = test.substr(found_position+1+a_isotope.length()); 01140 // Then we convert it into a string stream to ease the reading of both numbers 01141 istringstream fstr(test); 01142 fstr >> half_life >> branching_ratio; 01143 // We found it 01144 found_it = true; 01145 break; 01146 } 01147 // Read a new line 01148 fin.getline(line,line_max_size); 01149 } 01150 01151 delete[] line; 01152 // Close file 01153 fin.close(); 01154 01155 // Check if we found it or not 01156 if (found_it) 01157 { 01158 // Check rationality of values 01159 if (branching_ratio<=0. || branching_ratio>1.) 01160 { 01161 Cerr("***** oImageDimensionsAndQuantification::SetPETIsotope() -> Branching ratio (" << branching_ratio << ") is not in the ]0:1] range !" << endl); 01162 return 1; 01163 } 01164 // Verbose 01165 if (m_verbose>=2 && a_bed==m_nbBeds-1) 01166 { 01167 // If negative half-life, then we consider it is infinite (no half-life) 01168 if (half_life<=0.) 01169 Cout("oImageDimensionsAndQuantification::SetPETIsotope() -> Isotope " << a_isotope << " has infinite half life and " << branching_ratio << " branching ratio" << endl); 01170 else 01171 Cout("oImageDimensionsAndQuantification::SetPETIsotope() -> Isotope " << a_isotope << " has " << half_life << " seconds half life and " << branching_ratio << " branching ratio" << endl); 01172 } 01173 } 01174 else 01175 { 01176 // Throw error message 01177 Cerr("***** oImageDimensionsAndQuantification::SetPETIsotope() -> Did not find " << a_isotope << " isotope in the PET isotope data file, please add it !" << endl); 01178 return 1; 01179 } 01180 01181 // ------------------------------------------------------------------------ 01182 // Second step is applying branching ratio and decay wrt to frame durations 01183 // ------------------------------------------------------------------------ 01184 01185 // Branching ratio 01186 if (!m_ignoreBratCorrectionFlag) 01187 { 01188 // Apply correction 01189 for (int fr=0; fr<m_nbTimeFrames; fr++) 01190 { 01191 for(int g=0 ; g<m_nbRespGates*m_nbCardGates ; g++) 01192 m3p_quantificationFactors[a_bed][fr][g] /= branching_ratio; 01193 } 01194 // Verbose 01195 if (m_verbose>=2 && a_bed==m_nbBeds-1) Cout(" --> Correct for branching ratio" << endl); 01196 } 01197 // Verbose 01198 else if (m_verbose>=2 && a_bed==m_nbBeds-1) Cout(" --> Ignore branching ratio correction" << endl); 01199 01200 // Apply decay factors if half life is not negative (meaning no half-life) 01201 if (half_life>0.) 01202 { 01203 // We correct for it 01204 if (!m_ignoreDecaCorrectionFlag) 01205 { 01206 // Apply correction 01207 for (int fr=0; fr<m_nbTimeFrames; fr++) 01208 { 01209 FLTNB lambda = log(2.0)/half_life; 01210 for(int g=0 ; g<m_nbRespGates*m_nbCardGates ; g++) 01211 { 01212 // Time start decay correction 01213 m3p_quantificationFactors[a_bed][fr][g] *= exp(lambda*GetFrameTimeStartInSec(a_bed,fr)); 01214 // Frame duration decay correction 01215 m3p_quantificationFactors[a_bed][fr][g] *= lambda*GetFrameDurationInSec(a_bed,fr)/(1.0-exp(-lambda*(GetFrameDurationInSec(a_bed,fr)))); 01216 } 01217 } 01218 // Verbose 01219 if (m_verbose>=2 && a_bed==m_nbBeds-1) Cout(" --> Correct for half-life" << endl); 01220 } 01221 // Verbose 01222 else if (m_verbose>=2 && a_bed==m_nbBeds-1) Cout(" --> Ignore half-life correction" << endl); 01223 } 01224 01225 // End 01226 return 0; 01227 } 01228 01229 // ===================================================================== 01230 // --------------------------------------------------------------------- 01231 // --------------------------------------------------------------------- 01232 // ===================================================================== 01233 01234 int oImageDimensionsAndQuantification::InitDynamicData( string a_pathTo4DDataSplittingFile, 01235 int a_respMotionCorrectionFlag, 01236 int a_cardMotionCorrectionFlag, 01237 int a_invMotionCorrectionFlag, 01238 int a_nbRespGates, int a_nbCardGates ) 01239 { 01240 if (m_verbose>=5) Cout("oImageDimensionsAndQuantification::InitDynamicData()" << endl); 01241 mp_DynamicDataManager->SetVerbose(m_verbose); 01242 mp_DynamicDataManager->SetImageDimensionsAndQuantification(this); 01243 if (mp_DynamicDataManager->InitDynamicData( a_nbRespGates, a_nbCardGates, a_pathTo4DDataSplittingFile, 01244 a_respMotionCorrectionFlag, a_cardMotionCorrectionFlag, a_invMotionCorrectionFlag )) 01245 { 01246 Cerr("***** oImageDimensionsAndQuantification::InitDynamicData() -> A problem occured while initializing the dynamic data from dynamic data manager !" << endl); 01247 return 1; 01248 } 01249 return 0; 01250 } 01251 01252 // ===================================================================== 01253 // --------------------------------------------------------------------- 01254 // --------------------------------------------------------------------- 01255 // ===================================================================== 01256 01257 int oImageDimensionsAndQuantification::CheckDynamicParameters(int64_t a_nbEvents) 01258 { 01259 if (m_verbose>=5) Cout("oImageDimensionsAndQuantification::CheckDynamicParameters()" << endl); 01260 return mp_DynamicDataManager->CheckParameters(a_nbEvents); 01261 } 01262 01263 // ===================================================================== 01264 // --------------------------------------------------------------------- 01265 // --------------------------------------------------------------------- 01266 // =====================================================================