![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class iOptimizerMLEM 00004 00005 - separators: X 00006 - doxygen: X 00007 - default initialization: X 00008 - CASTOR_DEBUG: 00009 - CASTOR_VERBOSE: 00010 */ 00011 00018 #include "iOptimizerMLEM.hh" 00019 #include "sOutputManager.hh" 00020 00021 // ===================================================================== 00022 // --------------------------------------------------------------------- 00023 // --------------------------------------------------------------------- 00024 // ===================================================================== 00025 00026 iOptimizerMLEM::iOptimizerMLEM() : vOptimizer() 00027 { 00028 // --------------------------- 00029 // Mandatory member parameters 00030 // --------------------------- 00031 00032 // Initial value at 1 00033 m_initialValue = 1.; 00034 // Only one backward image for MLEM 00035 m_nbBackwardImages = 1; 00036 // MLEM does not accept penalties 00037 //m_penaltyEnergyFunctionDerivativesOrder = 0; 00038 // MLEM is compatible with listmode and histogram data 00039 m_listmodeCompatibility = true; 00040 m_histogramCompatibility = true; 00041 00042 // -------------------------- 00043 // Specific member parameters 00044 // -------------------------- 00045 00046 m_dataSpaceDenominatorThreshold = -1.; 00047 m_minimumImageUpdateFactor = -1.; 00048 m_maximumImageUpdateFactor = -1.; 00049 } 00050 00051 // ===================================================================== 00052 // --------------------------------------------------------------------- 00053 // --------------------------------------------------------------------- 00054 // ===================================================================== 00055 00056 iOptimizerMLEM::~iOptimizerMLEM() 00057 { 00058 } 00059 00060 // ===================================================================== 00061 // --------------------------------------------------------------------- 00062 // --------------------------------------------------------------------- 00063 // ===================================================================== 00064 00065 void iOptimizerMLEM::ShowHelpSpecific() 00066 { 00067 cout << "This optimizer is the standard MLEM (for Maximum Likelihood Expectation Maximization)." << endl; 00068 cout << "If subsets are used, it naturally becomes the OSEM optimizer." << endl; 00069 cout << "The following options can be used (in this particular order when provided as a list):" << endl; 00070 cout << " initial image value: to set the uniform voxel value for the initial image" << endl; 00071 cout << " denominator threshold: to set the threshold of the data space denominator under which the ratio is set to 1" << endl; 00072 cout << " minimum image update: to set the minimum of the image update factor under which it stays constant (0 or a negative value" << endl; 00073 cout << " means no minimum but still forbidding a 0 update)" << endl; 00074 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; 00075 cout << " no maximum)" << endl; 00076 } 00077 00078 // ===================================================================== 00079 // --------------------------------------------------------------------- 00080 // --------------------------------------------------------------------- 00081 // ===================================================================== 00082 00083 int iOptimizerMLEM::ReadConfigurationFile(const string& a_configurationFile) 00084 { 00085 string key_word = ""; 00086 // Read the initial image value option 00087 key_word = "initial image value"; 00088 if (ReadDataASCIIFile(a_configurationFile, key_word, &m_initialValue, 1, KEYWORD_MANDATORY)) 00089 { 00090 Cerr("***** iOptimizerMLEM::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl); 00091 return 1; 00092 } 00093 // Read the denominator threshold option 00094 key_word = "denominator threshold"; 00095 if (ReadDataASCIIFile(a_configurationFile, key_word, &m_dataSpaceDenominatorThreshold, 1, KEYWORD_MANDATORY)) 00096 { 00097 Cerr("***** iOptimizerMLEM::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl); 00098 return 1; 00099 } 00100 // Read the minimum image update option 00101 key_word = "minimum image update"; 00102 if (ReadDataASCIIFile(a_configurationFile, key_word, &m_minimumImageUpdateFactor, 1, KEYWORD_MANDATORY)) 00103 { 00104 Cerr("***** iOptimizerMLEM::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl); 00105 return 1; 00106 } 00107 // Read the maximum image update option 00108 key_word = "maximum image update"; 00109 if (ReadDataASCIIFile(a_configurationFile, key_word, &m_maximumImageUpdateFactor, 1, KEYWORD_MANDATORY)) 00110 { 00111 Cerr("***** iOptimizerMLEM::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl); 00112 return 1; 00113 } 00114 // Normal end 00115 return 0; 00116 } 00117 00118 // ===================================================================== 00119 // --------------------------------------------------------------------- 00120 // --------------------------------------------------------------------- 00121 // ===================================================================== 00122 00123 int iOptimizerMLEM::ReadOptionsList(const string& a_optionsList) 00124 { 00125 // There are 4 floating point variables as options 00126 FLTNB options[4]; 00127 00128 // Read them 00129 if (ReadStringOption(a_optionsList, options, 4, ",", "MLEM configuration")) 00130 { 00131 Cerr("***** iOptimizerMLEM::ReadAndCheckConfigurationFile() -> Failed to correctly read the list of options !" << endl); 00132 return 1; 00133 } 00134 // Affect options 00135 m_initialValue = options[0]; 00136 m_dataSpaceDenominatorThreshold = options[1]; 00137 m_minimumImageUpdateFactor = options[2]; 00138 m_maximumImageUpdateFactor = options[3]; 00139 00140 // Normal end 00141 return 0; 00142 } 00143 00144 // ===================================================================== 00145 // --------------------------------------------------------------------- 00146 // --------------------------------------------------------------------- 00147 // ===================================================================== 00148 00149 int iOptimizerMLEM::CheckSpecificParameters() 00150 { 00151 // Check that initial image value is strictly positive 00152 if (m_initialValue<=0.) 00153 { 00154 Cerr("***** iOptimizerMLEM->Initialize() -> Provided initial image value (" << m_initialValue << ") must be strictly positive !" << endl); 00155 return 1; 00156 } 00157 // Check that denominator threshold value is strictly positive 00158 if (m_dataSpaceDenominatorThreshold<=0.) 00159 { 00160 Cerr("***** iOptimizerMLEM->Initialize() -> Provided data space denominator threshold (" << m_dataSpaceDenominatorThreshold << ") must be strictly positive !" << endl); 00161 return 1; 00162 } 00163 // Check that maximum image update factor is higher than the minimum 00164 if (m_minimumImageUpdateFactor>0. && m_maximumImageUpdateFactor>0. && m_maximumImageUpdateFactor<m_minimumImageUpdateFactor) 00165 { 00166 Cerr("***** iOptimizerMLEM->Initialize() -> Provided minimum/maximum (" << m_minimumImageUpdateFactor << "/" << m_maximumImageUpdateFactor << " are inconsistent !" << endl); 00167 return 1; 00168 } 00169 // Normal end 00170 return 0; 00171 } 00172 00173 // ===================================================================== 00174 // --------------------------------------------------------------------- 00175 // --------------------------------------------------------------------- 00176 // ===================================================================== 00177 00178 int iOptimizerMLEM::InitializeSpecific() 00179 { 00180 // Verbose 00181 if (m_verbose>=1) 00182 { 00183 Cout("iOptimizerMLEM::Initialize() -> Use the MLEM optimizer" << endl); 00184 if (m_verbose>=2) 00185 { 00186 Cout(" --> Initial image value: " << m_initialValue << endl); 00187 Cout(" --> Data space denominator threshold: " << m_dataSpaceDenominatorThreshold << endl); 00188 if (m_minimumImageUpdateFactor>0.) Cout(" --> Minimum image update factor: " << m_minimumImageUpdateFactor << endl); 00189 else Cerr("!!!!! The minimum update value is not set, if using subsets, voxels could be trapped in 0 value causing some negative bias !" << endl); 00190 if (m_maximumImageUpdateFactor>0.) Cout(" --> Maximum image update factor: " << m_maximumImageUpdateFactor << endl); 00191 } 00192 } 00193 // Normal end 00194 return 0; 00195 } 00196 00197 // ===================================================================== 00198 // --------------------------------------------------------------------- 00199 // --------------------------------------------------------------------- 00200 // ===================================================================== 00201 00202 int iOptimizerMLEM::SensitivitySpecificOperations( FLTNB a_data, FLTNB a_forwardModel, FLTNB* ap_weight, 00203 FLTNB a_multiplicativeCorrections, FLTNB a_additiveCorrections, 00204 FLTNB a_quantificationFactor, oProjectionLine* ap_Line ) 00205 { 00206 // Line weight here is simply 1 00207 *ap_weight = 1.; 00208 // That's all 00209 return 0; 00210 } 00211 00212 // ===================================================================== 00213 // --------------------------------------------------------------------- 00214 // --------------------------------------------------------------------- 00215 // ===================================================================== 00216 00217 int iOptimizerMLEM::DataSpaceSpecificOperations( FLTNB a_data, FLTNB a_forwardModel, FLTNB* ap_backwardValues, 00218 FLTNB a_multiplicativeCorrections, FLTNB a_additiveCorrections, 00219 FLTNB a_quantificationFactor, oProjectionLine* ap_Line ) 00220 { 00221 // Truncate data to 0 if negative 00222 if (a_data<0.) a_data = 0.; 00223 // If the foward model is to close to zero, then ignore this data (set to 1 and not 0 because this line is taken into account in the sensitivity) 00224 if (a_forwardModel>m_dataSpaceDenominatorThreshold) *ap_backwardValues = a_data / a_forwardModel; 00225 else *ap_backwardValues = 1.; 00226 // That's all 00227 return 0; 00228 } 00229 00230 // ===================================================================== 00231 // --------------------------------------------------------------------- 00232 // --------------------------------------------------------------------- 00233 // ===================================================================== 00234 00235 int iOptimizerMLEM::ImageSpaceSpecificOperations( FLTNB a_currentImageValue, FLTNB* ap_newImageValue, 00236 FLTNB a_sensitivity, FLTNB* ap_correctionValues ) 00237 { 00238 // Compute image update factor 00239 FLTNB image_update_factor = *ap_correctionValues / a_sensitivity; 00240 // Apply minimum image update factor 00241 if ( m_minimumImageUpdateFactor > 0. && image_update_factor < m_minimumImageUpdateFactor ) image_update_factor = m_minimumImageUpdateFactor; 00242 // Apply maximum image update factor 00243 if ( m_maximumImageUpdateFactor > 0. && image_update_factor > m_maximumImageUpdateFactor ) image_update_factor = m_maximumImageUpdateFactor; 00244 // Update image 00245 *ap_newImageValue = a_currentImageValue * image_update_factor; 00246 // End 00247 return 0; 00248 } 00249 00250 // ===================================================================== 00251 // --------------------------------------------------------------------- 00252 // --------------------------------------------------------------------- 00253 // =====================================================================