103 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
104 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
105 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
106 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
107 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
108 Cerr(
"!!!!! iOptimizerOneStepLate::Destructor() -> Several times, the update factor was not positive and was set to 1. !!!!!" << endl);
109 Cerr(
"!!!!! For convenience, it was advertised only the first time. But this may be a sign that the penalty strength !!!!!" << endl);
110 Cerr(
"!!!!! is too high for this approximative algorithm. It may also be due to voxels close to FOV extremities. Be !!!!!" << endl);
111 Cerr(
"!!!!! sure to double check your images ! !!!!!" << endl);
112 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
113 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
114 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
115 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
116 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
127 cout <<
"This optimizer is the One-Step-Late algorithm from P. J. Green, IEEE TMI, Mar 1990, vol. 9, pp. 84-93." << endl;
128 cout <<
"Subsets can be used as for OSEM. It accepts penalty terms that have a derivative order of at least one." << endl;
129 cout <<
"Without penalty, it is stricly equivalent to the MLEM algorithm." << endl;
130 cout <<
"It is numerically implemented in the multiplicative form (as opposed to the gradient form)." << endl;
131 cout <<
"The following options can be used (in this particular order when provided as a list):" << endl;
132 cout <<
" initial image value: to set the uniform voxel value for the initial image" << endl;
133 cout <<
" denominator threshold: to set the threshold of the data space denominator under which the ratio is set to 1" << endl;
134 cout <<
" minimum image update: to set the minimum of the image update factor under which it stays constant (0 or a negative value" << endl;
135 cout <<
" means no minimum thus allowing a 0 update)" << endl;
136 cout <<
" maximum image update: to set the maximum of the image update factor over which it stays constant (0 or a negative value means" << endl;
137 cout <<
" no maximum)" << endl;
147 string key_word =
"";
149 key_word =
"initial image value";
152 Cerr(
"***** iOptimizerOneStepLate::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word <<
"' keyword !" << endl);
156 key_word =
"denominator threshold";
159 Cerr(
"***** iOptimizerOneStepLate::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word <<
"' keyword !" << endl);
163 key_word =
"minimum image update";
166 Cerr(
"***** iOptimizerOneStepLate::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word <<
"' keyword !" << endl);
170 key_word =
"maximum image update";
173 Cerr(
"***** iOptimizerOneStepLate::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word <<
"' keyword !" << endl);
188 const int nb_options = 4;
189 FLTNB options[nb_options];
191 if (
ReadStringOption(a_optionsList, options, nb_options,
",",
"OneStepLate configuration"))
193 Cerr(
"***** iOptimizerOneStepLate::ReadAndCheckConfigurationFile() -> Failed to correctly read the list of options !" << endl);
215 Cerr(
"***** iOptimizerOneStepLate->Initialize() -> Provided initial image value (" <<
m_initialValue <<
") must be strictly positive !" << endl);
221 Cerr(
"***** iOptimizerOneStepLate->Initialize() -> Provided data space denominator threshold (" <<
m_dataSpaceDenominatorThreshold <<
") must be strictly positive !" << endl);
233 Cerr(
"***** iOptimizerMLEM->CheckSpecificParameters() -> Cannot reconstruct list-mode transmission data !" << endl);
272 Cout(
"iOptimizerOneStepLate::Initialize() -> Use the OneStepLate optimizer" << endl);
278 else Cerr(
"!!!!! The minimum update value is not set, if using subsets, voxels could be trapped in 0 value causing some negative bias !" << endl);
301 if (
m_verbose>=1)
Cout(
"iOptimizerOneStepLate::PreImageUpdateSpecificStep() -> Compute penalty term" << endl);
306 Cerr(
"***** iOptimizerOneStepLate::PreImageUpdateSpecificStep() -> A problem occurred while computing the penalty pre-processing step !" << endl);
320 bool problem =
false;
324 #pragma omp parallel for private(v) schedule(guided) 330 th = omp_get_thread_num();
335 Cerr(
"***** iOptimizerOneStepLate::PreImageUpdateSpecificStep() -> A problem occurred while computing the penalty local pre-processing step for voxel " << v <<
" !" << endl);
344 Cerr(
"***** iOptimizerOneStepLate::PreImageUpdateSpecificStep() -> A problem occurred inside the multi-threaded loop, stop now !" << endl);
359 FLTNB a_multiplicativeCorrections,
FLTNB a_additiveCorrections,
FLTNB a_blankValue,
373 FLTNB a_multiplicativeCorrections,
FLTNB a_additiveCorrections,
FLTNB a_blankValue,
380 if (a_data<0.) a_data = 0.;
383 else *ap_backwardValues = 1.;
389 a_data -= a_additiveCorrections;
390 a_forwardModel -= a_additiveCorrections;
394 if (a_data<1. || a_forwardModel<1. || a_data>a_blankValue || a_forwardModel>a_blankValue)
396 *ap_backwardValues = 1.;
400 a_data = log(a_blankValue/a_data);
401 a_forwardModel = log(a_blankValue/a_forwardModel);
404 else *ap_backwardValues = 1.;
416 FLTNB a_sensitivity,
FLTNB* ap_correctionValues,
417 INTNB a_voxel,
int a_tbf,
int a_rbf,
int a_cbf )
420 FLTNB image_update_factor = 0.;
424 if (sensitivity_plus_penalty <= 0.)
427 image_update_factor = 1.;
430 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
431 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
432 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
433 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
434 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
435 Cerr(
"!!!!! iOptimizerOneStepLate::ImageSpaceSpecificOperations() -> The update factor was not positive ! It has been set to 1. !!!!!" << endl);
436 Cerr(
"!!!!! This may be a sign that the penalty strength is too high, but it may only be due to voxels close to FOV extremities. !!!!!" << endl);
437 Cerr(
"!!!!! For convenience, this warning will appear only for this update. Be careful about your results ! !!!!!" << endl);
438 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
439 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
440 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
441 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
442 Cerr(
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
451 image_update_factor = *ap_correctionValues / sensitivity_plus_penalty;
458 *ap_newImageValue = a_currentImageValue * image_update_factor;
FLTNB m_maximumImageUpdateFactor
int m_requiredPenaltyDerivativesOrder
int GetNbCardBasisFunctions()
Get the number of cardiac basis functions.
Declaration of class iOptimizerOneStepLate.
void ShowHelpSpecific()
A function used to show help about the child optimizer.
bool m_listmodeCompatibility
int SensitivitySpecificOperations(FLTNB a_data, FLTNB a_forwardModel, FLTNB *ap_weight, FLTNB a_multiplicativeCorrections, FLTNB a_additiveCorrections, FLTNB a_blankValue, FLTNB a_quantificationFactor, oProjectionLine *ap_Line)
This function compute the weight associated to the provided event (for sensitivity computation) ...
oImageDimensionsAndQuantification * mp_ImageDimensionsAndQuantification
FLTNB **** m4p_firstDerivativePenaltyImage
virtual int LocalPreProcessingStep(int a_tbf, int a_rbf, int a_cbf, INTNB a_voxel, int a_th)
A public function computing a local pre-processing step for the penalty.
int GetNbTimeBasisFunctions()
Get the number of time basis functions.
int ImageSpaceSpecificOperations(FLTNB a_currentImageValue, FLTNB *ap_newImageValue, FLTNB a_sensitivity, FLTNB *ap_correctionValues, INTNB a_voxel, int a_tbf=-1, int a_rbf=-1, int a_cbf=-1)
This function perform the image update step specific to the optimizer.
int ReadOptionsList(const string &a_optionsList)
A function used to read options from a list of options.
bool m_emissionCompatibility
int ReadConfigurationFile(const string &a_configurationFile)
A function used to read options from a configuration file.
bool m_displayWarningFlag
virtual int GlobalPreProcessingStep()
A public function computing a global pre-processing step for the penalty.
int ReadDataASCIIFile(const string &a_file, const string &a_keyword, T *ap_return, int a_nbElts, bool a_mandatoryFlag)
Look for "a_nbElts" elts in the "a_file" file matching the "a_keyword" string passed as parameter a...
FLTNB m_minimumImageUpdateFactor
~iOptimizerOneStepLate()
The destructor of iOptimizerOneStepLate.
bool m_histogramCompatibility
#define KEYWORD_MANDATORY
#define SPEC_TRANSMISSION
bool m_transmissionCompatibility
This class is designed to generically described any iterative optimizer.
This class is designed to manage and store system matrix elements associated to a vEvent...
int DataSpaceSpecificOperations(FLTNB a_data, FLTNB a_forwardModel, FLTNB *ap_backwardValues, FLTNB a_multiplicativeCorrections, FLTNB a_additiveCorrections, FLTNB a_blankValue, FLTNB a_quantificationFactor, oProjectionLine *ap_Line)
This function performs the data space operations specific to the optimizer (computes the values to be...
Declaration of class sOutputManager.
oImageSpace * mp_ImageSpace
int GetNbThreadsForImageComputation()
Get the number of threads used for image operations.
FLTNB m_dataSpaceDenominatorThreshold
INTNB GetNbVoxXYZ()
Get the total number of voxels.
iOptimizerOneStepLate()
The constructor of iOptimizerOneStepLate.
int InitializeSpecific()
This function is used to initialize specific stuff to the child optimizer.
int ReadStringOption(const string &a_input, T *ap_return, int a_nbElts, const string &sep, const string &a_option)
Parse the 'a_input' string corresponding to the 'a_option' into 'a_nbElts' elements, using the 'sep' separator. The results are returned in the templated 'ap_return' dynamic templated array. Call "ConvertFromString()" to perform the correct conversion depending on the type of the data to convert.
int GetNbRespBasisFunctions()
Get the number of respiratory basis functions.
virtual FLTNB ComputeFirstDerivative(int a_tbf, int a_rbf, int a_cbf, INTNB a_voxel, int a_th)=0
A public function computing the derivative of the penalty.
int PreImageUpdateSpecificStep()
A private function used to compute the penalty term of the OneStepLate algorithm. ...
int CheckSpecificParameters()
A private function used to check the parameters settings specific to the child optimizer.