![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class oOptimizerManager 00004 00005 - separators: X 00006 - doxygen: X 00007 - default initialization: X 00008 - CASTOR_DEBUG: 00009 - CASTOR_VERBOSE: 00010 */ 00011 00018 #include "oOptimizerManager.hh" 00019 #include "sOutputManager.hh" 00020 #include "sAddonManager.hh" 00021 00022 // ===================================================================== 00023 // --------------------------------------------------------------------- 00024 // --------------------------------------------------------------------- 00025 // ===================================================================== 00026 00027 oOptimizerManager::oOptimizerManager() 00028 { 00029 // Image dimensions 00030 mp_ImageDimensionsAndQuantification = NULL; 00031 // Data mode 00032 m_dataMode = MODE_UNKNOWN; 00033 // Data type 00034 m_dataType = TYPE_UNKNOWN; 00035 // Number of TOF bins 00036 m_nbTOFBins = 0; 00037 // Optimizer and penalty options 00038 m_optionsOptimizer = ""; 00039 m_optionsPenalty = ""; 00040 // Optimizer and penalty 00041 mp_Optimizer = NULL; 00042 mp_Penalty = NULL; 00043 // Verbosity 00044 m_verbose = 0; 00045 // Optimizer FOM computation, and image update stat flags 00046 m_optimizerFOMFlag = false; 00047 m_optimizerImageStatFlag = false; 00048 } 00049 00050 // ===================================================================== 00051 // --------------------------------------------------------------------- 00052 // --------------------------------------------------------------------- 00053 // ===================================================================== 00054 00055 oOptimizerManager::~oOptimizerManager() 00056 { 00057 } 00058 00059 // ===================================================================== 00060 // --------------------------------------------------------------------- 00061 // --------------------------------------------------------------------- 00062 // ===================================================================== 00063 00064 int oOptimizerManager::CheckParameters() 00065 { 00066 // Check image dimensions 00067 if (mp_ImageDimensionsAndQuantification==NULL) 00068 { 00069 Cerr("***** oOptimizerManager::CheckParameters() -> No image dimensions provided !" << endl); 00070 return 1; 00071 } 00072 // Check data mode 00073 if (m_dataMode==MODE_UNKNOWN || (m_dataMode!=MODE_LIST && m_dataMode!=MODE_HISTOGRAM)) 00074 { 00075 Cerr("***** oOptimizerManager::CheckParameters() -> No or meaningless data mode provided !" << endl); 00076 return 1; 00077 } 00078 // Check data type 00079 if (m_dataType==TYPE_UNKNOWN || (m_dataType!=TYPE_PET && m_dataType!=TYPE_SPECT && m_dataType!=TYPE_TRANSMISSION)) 00080 { 00081 Cerr("***** oOptimizerManager::CheckParameters() -> No or meaningless data type provided !" << endl); 00082 return 1; 00083 } 00084 // Check optimizer options 00085 if (m_optionsOptimizer=="") 00086 { 00087 Cerr("***** oOptimizerManager::CheckParameters() -> No optimizer options provided !" << endl); 00088 return 1; 00089 } 00090 // Check verbosity 00091 if (m_verbose<0) 00092 { 00093 Cerr("***** oOptimizerManager::CheckParameters() -> Wrong verbosity level provided !" << endl); 00094 return 1; 00095 } 00096 // Normal end 00097 return 0; 00098 } 00099 00100 // ===================================================================== 00101 // --------------------------------------------------------------------- 00102 // --------------------------------------------------------------------- 00103 // ===================================================================== 00104 00105 int oOptimizerManager::Initialize() 00106 { 00107 // Verbose 00108 if (m_verbose>=1) Cout("oOptimizerManager::Initialize() -> Initialize optimizer and penalty" << endl); 00109 00110 // Parse projector options and initialize them 00111 if (ParseOptionsAndInitializeOptimizerAndPenalty()) 00112 { 00113 Cerr("***** oOptimizerManager::Initialize() -> A problem occured while parsing optimizer options and initializing it !" << endl); 00114 return 1; 00115 } 00116 00117 // Normal end 00118 return 0; 00119 } 00120 00121 // ===================================================================== 00122 // --------------------------------------------------------------------- 00123 // --------------------------------------------------------------------- 00124 // ===================================================================== 00125 00126 int oOptimizerManager::ParseOptionsAndInitializeOptimizerAndPenalty() 00127 { 00128 // --------------------------------------------------------------------------------------------------- 00129 // Manage optimizer options 00130 // --------------------------------------------------------------------------------------------------- 00131 00132 string name_optimizer = ""; 00133 string list_options_optimizer = ""; 00134 string file_options_optimizer = ""; 00135 00136 // This is for the automatic initialization of the optimizer and penalty 00137 typedef vOptimizer *(*maker_optimizer) (); 00138 00139 // ______________________________________________________________________________ 00140 // Get the optimizer name in the options and isolate the real optimizer's options 00141 00142 // First check emptyness of the options 00143 if (m_optionsOptimizer=="") 00144 { 00145 Cerr("***** oOptimizerManager::ParseOptionsAndInitializeOptimizerAndPenalty() -> No optimizer provided !" << endl); 00146 return 1; 00147 } 00148 00149 // Search for a colon ":", this indicates that a configuration file is provided after the optimizer name 00150 size_t colon = m_optionsOptimizer.find_first_of(":"); 00151 size_t comma = m_optionsOptimizer.find_first_of(","); 00152 00153 // Case 1: we have a colon 00154 if (colon!=string::npos) 00155 { 00156 // Get the optimizer name before the colon 00157 name_optimizer = m_optionsOptimizer.substr(0,colon); 00158 // Get the configuration file after the colon 00159 file_options_optimizer = m_optionsOptimizer.substr(colon+1); 00160 // List of options is empty 00161 list_options_optimizer = ""; 00162 } 00163 // Case 2: we have a comma 00164 else if (comma!=string::npos) 00165 { 00166 // Get the optimizer name before the first comma 00167 name_optimizer = m_optionsOptimizer.substr(0,comma); 00168 // Get the list of options after the first comma 00169 list_options_optimizer = m_optionsOptimizer.substr(comma+1); 00170 // Configuration file is empty 00171 file_options_optimizer = ""; 00172 } 00173 // Case 3: no colon and no comma (a single optimizer name) 00174 else 00175 { 00176 // Get the optimizer name 00177 name_optimizer = m_optionsOptimizer; 00178 // List of options is empty 00179 list_options_optimizer = ""; 00180 // Build the default configuration file 00181 sOutputManager* p_output_manager = sOutputManager::GetInstance(); 00182 file_options_optimizer = p_output_manager->GetPathToConfigDir() + "/optimizer/" + name_optimizer + ".conf"; 00183 } 00184 00185 // ______________________________________________________________________________ 00186 // Construct and initialize the optimizer 00187 00188 // Get optimizer's listfrom addon manager 00189 std::map <string,maker_optimizer> list_optimizer = sAddonManager::GetInstance()->mp_listOfOptimizers; 00190 // Create the optimizer 00191 if (list_optimizer[name_optimizer]) mp_Optimizer = list_optimizer[name_optimizer](); 00192 else 00193 { 00194 Cerr("***** oOptimizerManager::ParseOptionsAndInitializeOptimizerAndPenalty() -> Optimizer '" << name_optimizer << "' does not exist !" << endl); 00195 return 1; 00196 } 00197 // Set parameters 00198 mp_Optimizer->SetImageDimensionsAndQuantification(mp_ImageDimensionsAndQuantification); 00199 mp_Optimizer->SetDataMode(m_dataMode); 00200 mp_Optimizer->SetDataType(m_dataType); 00201 mp_Optimizer->SetNbTOFBins(m_nbTOFBins); 00202 mp_Optimizer->SetFOMFlag(m_optimizerFOMFlag); 00203 mp_Optimizer->SetImageStatFlag(m_optimizerImageStatFlag); 00204 mp_Optimizer->SetVerbose(m_verbose); 00205 // Provide configuration file if any (child specific function) 00206 if (file_options_optimizer!="" && mp_Optimizer->ReadConfigurationFile(file_options_optimizer)) 00207 { 00208 Cerr("***** oOptimizerManager::ParseOptionsAndInitializeOptimizerAndPenalty() -> A problem occured while reading and checking optimizer's configuration file !" << endl); 00209 return 1; 00210 } 00211 // Provide options if any (child specific function) 00212 if (list_options_optimizer!="" && mp_Optimizer->ReadOptionsList(list_options_optimizer)) 00213 { 00214 Cerr("***** oOptimizerManager::ParseOptionsAndInitializeOptimizerAndPenalty() -> A problem occured while parsing and reading optimizer's options !" << endl); 00215 return 1; 00216 } 00217 // Check parameters (mother generic function that will call the child specific function at the end) 00218 if (mp_Optimizer->CheckParameters()) 00219 { 00220 Cerr("***** oOptimizerManager::ParseOptionsAndInitializeOptimizerAndPenalty() -> A problem occured while checking optimizer parameters !" << endl); 00221 return 1; 00222 } 00223 // Initialize the optimizer (mother generic function that will call the child specific function at the end) 00224 if (mp_Optimizer->Initialize()) 00225 { 00226 Cerr("***** oOptimizerManager::ParseOptionsAndInitializeOptimizerAndPenalty() -> A problem occured while initializing the optimizer !" << endl); 00227 return 1; 00228 } 00229 00230 // --------------------------------------------------------------------------------------------------- 00231 // Manage penalty options (not yet implemented, a lot of work to do) 00232 // --------------------------------------------------------------------------------------------------- 00233 /* Do not remove 00234 string name_penalty = ""; 00235 string list_options_penalty = ""; 00236 string file_options_penalty = ""; 00237 00238 // This is for the automatic initialization of the penalty 00239 typedef vPenalty *(*maker_penalty) (); 00240 00241 // ______________________________________________________________________________ 00242 // Get the penalty name in the options and isolate the real penalty's options 00243 00244 // If no penalty, then return now 00245 if (m_optionsPenalty=="") 00246 { 00247 return 0; 00248 } 00249 // Otherwise, check if the optimizer accepts penalties 00250 else 00251 { 00252 if (!mp_Optimizer->GetAcceptPenalty()) 00253 { 00254 Cerr("***** oOptimizerManager::ParseOptionsAndInitializeOptimizerAndPenalty() -> Penalty provided while the selected optimizer does not accept penalties !" << endl); 00255 Cerr(" Remove penalty or change for another optimizer that accepts penalties." << endl); 00256 return 1; 00257 } 00258 } 00259 00260 // Search for a colon ":", this indicates that a configuration file is provided after the penalty name 00261 colon = m_optionsPenalty.find_first_of(":"); 00262 comma = m_optionsPenalty.find_first_of(","); 00263 00264 // Case 1: we have a colon 00265 if (colon!=string::npos) 00266 { 00267 // Get the penalty name before the colon 00268 name_penalty = m_optionsPenalty.substr(0,colon); 00269 // Get the configuration file after the colon 00270 file_options_penalty = m_optionsPenalty.substr(colon+1); 00271 // List of options is empty 00272 list_options_penalty = ""; 00273 } 00274 // Case 2: we have a comma 00275 else if (comma!=string::npos) 00276 { 00277 // Get the penalty name before the first comma 00278 name_penalty = m_optionsPenalty.substr(0,comma); 00279 // Get the list of options after the first comma 00280 list_options_penalty = m_optionsPenalty.substr(comma+1); 00281 // Configuration file is empty 00282 file_options_penalty = ""; 00283 } 00284 // Case 3: no colon and no comma (a single penalty name) 00285 else 00286 { 00287 // Get the penalty name 00288 name_penalty = m_optionsPenalty; 00289 // List of options is empty 00290 list_options_penalty = ""; 00291 // Build the default configuration file 00292 sOutputManager* p_output_manager = sOutputManager::GetInstance(); 00293 file_options_penalty = p_output_manager->GetPathToConfigDir() + "/penalty/" + name_optimizer + ".conf"; 00294 } 00295 00296 // ______________________________________________________________________________ 00297 // Construct and initialize the penalty 00298 00299 // Get penalty's list from addon manager 00300 std::map <string,maker_penalty> list_penalty = sAddonManager::GetInstance()->mp_listOfPenalties; 00301 // Create the penalty 00302 if (list_penalty[name_penalty]) mp_Penalty = list_penalty[name_penalty](); 00303 else 00304 { 00305 Cerr("***** oOptimizerManager::ParseOptionsAndInitializeOptimizerAndPenalty() -> Penalty '" << name_penalty << "' does not exist !" << endl); 00306 return 1; 00307 } 00308 // Set parameters 00309 mp_Penalty->SetImageDimensionsAndQuantification(mp_ImageDimensionsAndQuantification); 00310 mp_Penalty->SetVerbose(m_verbose); 00311 // Check parameters 00312 if (mp_Penalty->CheckParameters()) 00313 { 00314 Cerr("***** oOptimizerManager::ParseOptionsAndInitializeOptimizerAndPenalty() -> A problem occured while checking penalty parameters !" << endl); 00315 return 1; 00316 } 00317 // Provide configuration file if any 00318 if (file_options_penalty!="" && mp_Penalty->ReadAndCheckConfigurationFile(file_options_penalty)) 00319 { 00320 Cerr("***** oOptimizerManager::ParseOptionsAndInitializeOptimizerAndPenalty() -> A problem occured while reading and checking penalty's configuration file !" << endl); 00321 return 1; 00322 } 00323 // Provide options if any 00324 if (list_options_penalty!="" && mp_Penalty->ReadAndCheckOptionsList(list_options_penalty)) 00325 { 00326 Cerr("***** oOptimizerManager::ParseOptionsAndInitializeOptimizerAndPenalty() -> A problem occured while parsing and reading penalty's options !" << endl); 00327 return 1; 00328 } 00329 // Initialize the penalty 00330 if (mp_Penalty->Initialize()) 00331 { 00332 Cerr("***** oOptimizerManager::ParseOptionsAndInitializeOptimizerAndPenalty() -> A problem occured while initializing the penalty !" << endl); 00333 return 1; 00334 } 00335 // Check that the derivatives order of the penalty is compatible with the optimizer 00336 if (mp_Penalty->GetPenaltyEnergyFunctionDerivativesOrder() != mp_Optimizer->GetPenaltyEnergyFunctionDerivativesOrder()) 00337 { 00338 Cerr("***** oOptimizerManager::ParseOptionsAndInitializeOptimizerAndPenalty() -> Derivatives order allowed by chosen penalty is not compatible with the order accepted by the optimizer !" << endl); 00339 return 1; 00340 } 00341 */ 00342 00343 // Normal end 00344 return 0; 00345 } 00346 00347 // ===================================================================== 00348 // --------------------------------------------------------------------- 00349 // --------------------------------------------------------------------- 00350 // ===================================================================== 00351 00352 int oOptimizerManager::PreDataUpdateStep(int a_iteration, int a_nbIterations, int a_subset, int a_nbSubsets) 00353 { 00354 // Simply call the homonymous function from the vOptimizer 00355 if (mp_Optimizer->PreDataUpdateStep(a_iteration, a_nbIterations, a_subset, a_nbSubsets)) 00356 { 00357 Cerr("***** oOptimizerManager::PreDataUpdateStep() -> A problem occured while applying the pre-data-update step to the optimizer !" << endl); 00358 return 1; 00359 } 00360 // Normal end 00361 return 0; 00362 } 00363 00364 // ===================================================================== 00365 // --------------------------------------------------------------------- 00366 // --------------------------------------------------------------------- 00367 // ===================================================================== 00368 00369 int oOptimizerManager::PostDataUpdateStep(int a_iteration, int a_nbIterations, int a_subset, int a_nbSubsets) 00370 { 00371 // Simply call the homonymous function from the vOptimizer 00372 if (mp_Optimizer->PostDataUpdateStep(a_iteration, a_nbIterations, a_subset, a_nbSubsets)) 00373 { 00374 Cerr("***** oOptimizerManager::PostDataUpdateStep() -> A problem occured while applying the post-data-update step to the optimizer !" << endl); 00375 return 1; 00376 } 00377 // Normal end 00378 return 0; 00379 } 00380 00381 // ===================================================================== 00382 // --------------------------------------------------------------------- 00383 // --------------------------------------------------------------------- 00384 // ===================================================================== 00385 00386 int oOptimizerManager::DataUpdateStep( oProjectionLine* ap_Line, oImageSpace* ap_Image, vEvent* ap_Event, 00387 int a_bed, int a_timeFrame, int a_respGate, int a_cardGate, 00388 int a_iteration, int a_thread) 00389 { 00390 // --------------------------------------------------------------------------------- 00391 // Deal with all multiplicative correction factors to be included in the projections 00392 // --------------------------------------------------------------------------------- 00393 00394 // Compute the global multiplicative correction factor 00395 FLTNB multiplicative_correction = ap_Event->GetMultiplicativeCorrections() * mp_ImageDimensionsAndQuantification->GetQuantificationFactor(a_bed,a_timeFrame, a_respGate, a_cardGate); 00396 00397 // Do nothing if the multiplicative correction factor is null or negative 00398 if (multiplicative_correction<=0.) return 0; 00399 00400 // Set the multiplicative correction to the oProjectionLine so that it is part of the system matrix and automatically applied 00401 ap_Line->SetMultiplicativeCorrection(multiplicative_correction); 00402 00403 // Set the current attenuation image for SPECT attenuation correction 00404 if (m_dataType==TYPE_SPECT && ap_Image->m4p_attenuation!=NULL) mp_Optimizer->SetAttenuationImage(ap_Image->m4p_attenuation[a_timeFrame][a_respGate][a_cardGate], a_thread); 00405 00406 // -------------------------------------------------------------------------------- 00407 // Decompose the data update step into 4 steps mandatory steps and 3 optional steps 00408 // -------------------------------------------------------------------------------- 00409 00410 // Mandatory 1: Compute model (forward-projection + additive background) 00411 if (mp_Optimizer->DataStep1ForwardProjectModel( ap_Line, ap_Image, ap_Event, a_bed, a_timeFrame, a_respGate, a_cardGate, a_thread )) 00412 { 00413 Cerr("***** oOptimizerManager::DataUpdateStep() -> An error occured while forward projecting !" << endl); 00414 return 1; 00415 } 00416 00417 // Optional 1: Do what is done in the child optimizer 00418 if (mp_Optimizer->DataStep2Optional( ap_Line, ap_Image, ap_Event, a_bed, a_timeFrame, a_respGate, a_cardGate, a_iteration, a_thread )) 00419 { 00420 Cerr("***** oOptimizerManager::DataUpdateStep() -> An error occured while performing optional step 1 !" << endl); 00421 return 1; 00422 } 00423 00424 // Mandatory 2: Compute sensitivity in histogram mode 00425 if (ap_Event->GetDataMode()==MODE_HISTOGRAM && mp_Optimizer->DataStep3BackwardProjectSensitivity( ap_Line, ap_Image, ap_Event, a_bed, a_timeFrame, a_respGate, a_cardGate, a_thread )) 00426 { 00427 Cerr("***** oOptimizerManager::DataUpdateStep() -> An error occured while backward projecting the sensitivity !" << endl); 00428 return 1; 00429 } 00430 00431 // Optional 2: Do what is done in the child optimizer 00432 if (mp_Optimizer->DataStep4Optional( ap_Line, ap_Image, ap_Event, a_bed, a_timeFrame, a_respGate, a_cardGate, a_iteration, a_thread )) 00433 { 00434 Cerr("***** oOptimizerManager::DataUpdateStep() -> An error occured while performing optional step 2 !" << endl); 00435 return 1; 00436 } 00437 00438 // Mandatory 3: Compute the correction terms 00439 if (mp_Optimizer->DataStep5ComputeCorrections( ap_Line, ap_Event, a_bed, a_timeFrame, a_respGate, a_cardGate, a_thread )) 00440 { 00441 Cerr("***** oOptimizerManager::DataUpdateStep() -> An error occured while computing correction terms !" << endl); 00442 return 1; 00443 } 00444 00445 // Optional 3: Do what is done in the child optimizer 00446 if (mp_Optimizer->DataStep6Optional( ap_Line, ap_Image, ap_Event, a_bed, a_timeFrame, a_respGate, a_cardGate, a_iteration, a_thread )) 00447 { 00448 Cerr("***** oOptimizerManager::DataUpdateStep() -> An error occured while performing optional step 3 !" << endl); 00449 return 1; 00450 } 00451 00452 // Mandatory 4: Backproject correction terms 00453 if (mp_Optimizer->DataStep7BackwardProjectCorrections( ap_Line, ap_Image, ap_Event, a_bed, a_timeFrame, a_respGate, a_cardGate, a_thread )) 00454 { 00455 Cerr("***** oOptimizerManager::DataUpdateStep() -> An error occured while backward projecting the correction !" << endl); 00456 return 1; 00457 } 00458 00459 // Compute FOM is asked for 00460 if (m_optimizerFOMFlag && mp_Optimizer->DataStep8ComputeFOM( ap_Line, ap_Event, a_timeFrame, a_respGate, a_cardGate, a_thread )) 00461 { 00462 Cerr("***** oOptimizerManager::DataUpdateStep() -> An error occured while computing FOMs !" << endl); 00463 return 1; 00464 } 00465 00466 // End 00467 return 0; 00468 } 00469 00470 // ===================================================================== 00471 // --------------------------------------------------------------------- 00472 // --------------------------------------------------------------------- 00473 // ===================================================================== 00474 00475 int oOptimizerManager::ImageUpdateStep(oImageSpace* ap_Image, int a_iteration, int a_nbSubsets) 00476 { 00477 // Verbose 00478 if (m_verbose>=2) Cout("oOptimizerManager::ImageUpdateStep() -> Proceed to image update" << endl); 00479 00480 // Update visited voxels 00481 if (mp_Optimizer->UpdateVisitedVoxels( ap_Image )) 00482 { 00483 Cerr("***** oOptimizerManager::ImageUpdateStep() -> Problem while updating visited voxels !" << endl); 00484 return 1; 00485 } 00486 00487 // Manage penalty computation 00488 if (mp_Penalty!=NULL) 00489 { 00490 // Loops on frames, etc ... 00491 } 00492 00493 // Image update step 00494 if (mp_Optimizer->ImageUpdateStep( ap_Image, a_iteration, a_nbSubsets )) 00495 { 00496 Cerr("***** oOptimizerManager::ImageUpdateStep() -> Problem while updating image space !" << endl); 00497 return 1; 00498 } 00499 00500 // End 00501 return 1; 00502 } 00503 00504 // ===================================================================== 00505 // --------------------------------------------------------------------- 00506 // --------------------------------------------------------------------- 00507 // =====================================================================