![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class iOptimizerOriginalAML 00004 00005 - separators: X 00006 - doxygen: X 00007 - default initialization: X 00008 - CASTOR_DEBUG: 00009 - CASTOR_VERBOSE: 00010 */ 00011 00018 #include "iOptimizerOriginalAML.hh" 00019 #include "sOutputManager.hh" 00020 00021 // ===================================================================== 00022 // --------------------------------------------------------------------- 00023 // --------------------------------------------------------------------- 00024 // ===================================================================== 00025 00026 iOptimizerOriginalAML::iOptimizerOriginalAML() : vOptimizer() 00027 { 00028 // --------------------------- 00029 // Mandatory member parameters 00030 // --------------------------- 00031 00032 // Initial value at 1 00033 m_initialValue = 1.; 00034 // Only one backward image for AML 00035 m_nbBackwardImages = 1; 00036 // AML does not accept penalties 00037 //m_penaltyEnergyFunctionDerivativesOrder = 0; 00038 // AML is only compatible with histogram data 00039 m_listmodeCompatibility = false; 00040 m_histogramCompatibility = true; 00041 00042 // -------------------------- 00043 // Specific member parameters 00044 // -------------------------- 00045 00046 // Threshold applied to the denominator of the data space operation to avoid to high ratios 00047 m_dataSpaceDenominatorThreshold = -1.; 00048 // This is the AML bound parameter 00049 m_bound = 1.; 00050 } 00051 00052 // ===================================================================== 00053 // --------------------------------------------------------------------- 00054 // --------------------------------------------------------------------- 00055 // ===================================================================== 00056 00057 iOptimizerOriginalAML::~iOptimizerOriginalAML() 00058 { 00059 } 00060 00061 // ===================================================================== 00062 // --------------------------------------------------------------------- 00063 // --------------------------------------------------------------------- 00064 // ===================================================================== 00065 00066 void iOptimizerOriginalAML::ShowHelpSpecific() 00067 { 00068 cout << "This optimizer is the AML algorithm derived from the AB-EMML of C. Byrne, Inverse Problems, 1998, vol. 14, pp. 1455-67." << endl; 00069 cout << "The bound B is taken to infinity, so only the bound A can be parameterized." << endl; 00070 cout << "This bound must be quantitative (same unit as the reconstructed image)." << endl; 00071 cout << "It is provided as a single value and thus assuming a uniform bound." << endl; 00072 cout << "Subsets can be used." << endl; 00073 cout << "With a negative or null bound, this algorithm implements equation 6 of A. Rahmim et al, Phys. Med. Biol., 2012, vol. 57, pp. 733-55." << endl; 00074 cout << "If a positive bound is provided, then we suppose that the bound A is taken to minus infinity. In that case, this algorithm implements" << endl; 00075 cout << "equation 22 of K. Van Slambrouck et al, IEEE TMI, Jan 2015, vol. 34, pp. 126-136." << endl; 00076 cout << "The following options can be used (in this particular order when provided as a list):" << endl; 00077 cout << " initial image value: to set the uniform voxel value for the initial image" << endl; 00078 cout << " denominator threshold: to set the threshold of the data space denominator under which the ratio is set to 1" << endl; 00079 cout << " bound: to set the bound parameter that shift the Poisson law (quantitative, negative or null for standard AML and positive for infinite AML)." << endl; 00080 } 00081 00082 // ===================================================================== 00083 // --------------------------------------------------------------------- 00084 // --------------------------------------------------------------------- 00085 // ===================================================================== 00086 00087 int iOptimizerOriginalAML::ReadConfigurationFile(const string& a_configurationFile) 00088 { 00089 string key_word = ""; 00090 // Read the initial image value option 00091 key_word = "initial image value"; 00092 if (ReadDataASCIIFile(a_configurationFile, key_word, &m_initialValue, 1, KEYWORD_MANDATORY)) 00093 { 00094 Cerr("***** iOptimizerOriginalAML::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl); 00095 return 1; 00096 } 00097 // Read the denominator threshold option 00098 key_word = "denominator threshold"; 00099 if (ReadDataASCIIFile(a_configurationFile, key_word, &m_dataSpaceDenominatorThreshold, 1, KEYWORD_MANDATORY)) 00100 { 00101 Cerr("***** iOptimizerOriginalAML::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl); 00102 return 1; 00103 } 00104 // Read the bound value 00105 key_word = "bound"; 00106 if (ReadDataASCIIFile(a_configurationFile, key_word, &m_bound, 1, KEYWORD_MANDATORY)) 00107 { 00108 Cerr("***** iOptimizerOriginalAML::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl); 00109 return 1; 00110 } 00111 // Normal end 00112 return 0; 00113 } 00114 00115 // ===================================================================== 00116 // --------------------------------------------------------------------- 00117 // --------------------------------------------------------------------- 00118 // ===================================================================== 00119 00120 int iOptimizerOriginalAML::ReadOptionsList(const string& a_optionsList) 00121 { 00122 // There are 3 floating point variables as options 00123 FLTNB options[3]; 00124 // Read them 00125 if (ReadStringOption(a_optionsList, options, 3, ",", "AML configuration")) 00126 { 00127 Cerr("***** iOptimizerOriginalAML::ReadAndCheckConfigurationFile() -> Failed to correctly read the list of options !" << endl); 00128 return 1; 00129 } 00130 // Affect options 00131 m_initialValue = options[0]; 00132 m_dataSpaceDenominatorThreshold = options[1]; 00133 m_bound = options[2]; 00134 // Normal end 00135 return 0; 00136 } 00137 00138 // ===================================================================== 00139 // --------------------------------------------------------------------- 00140 // --------------------------------------------------------------------- 00141 // ===================================================================== 00142 00143 int iOptimizerOriginalAML::CheckSpecificParameters() 00144 { 00145 // Check that denominator threshold value is strictly positive 00146 if (m_dataSpaceDenominatorThreshold<=0.) 00147 { 00148 Cerr("***** iOptimizerOriginalAML->Initialize() -> Provided data space denominator threshold (" << m_dataSpaceDenominatorThreshold << ") must be strictly positive !" << endl); 00149 return 1; 00150 } 00151 // If a negative or null bound is provided (standard AML), then check that the initial image value is above the provided bound 00152 if (m_bound<=0. && m_initialValue<m_bound) 00153 { 00154 Cerr("***** iOptimizerOriginalAML::Initialize() -> The initial image value (" << m_initialValue << ") must be higher than or equal to the provided bound value (" 00155 << m_bound << ") !" << endl); 00156 return 1; 00157 } 00158 // Normal end 00159 return 0; 00160 } 00161 00162 // ===================================================================== 00163 // --------------------------------------------------------------------- 00164 // --------------------------------------------------------------------- 00165 // ===================================================================== 00166 00167 int iOptimizerOriginalAML::InitializeSpecific() 00168 { 00169 // Verbose 00170 if (m_verbose>=1) 00171 { 00172 Cout("iOptimizerOriginalAML::Initialize() -> Use the AML optimizer" << endl); 00173 if (m_verbose>=2) 00174 { 00175 Cout(" --> Initial image value: " << m_initialValue << endl); 00176 Cout(" --> Data space denominator threshold: " << m_dataSpaceDenominatorThreshold << endl); 00177 if (m_bound<=0.) Cout(" --> Bound value: " << m_bound << endl); 00178 else Cout(" --> Bound value taken to minus infinity" << endl); 00179 } 00180 } 00181 // Normal end 00182 return 0; 00183 } 00184 00185 // ===================================================================== 00186 // --------------------------------------------------------------------- 00187 // --------------------------------------------------------------------- 00188 // ===================================================================== 00189 00190 int iOptimizerOriginalAML::SensitivitySpecificOperations( FLTNB a_data, FLTNB a_forwardModel, FLTNB* ap_weight, 00191 FLTNB a_multiplicativeCorrections, FLTNB a_additiveCorrections, 00192 FLTNB a_quantificationFactor, oProjectionLine* ap_Line ) 00193 { 00194 // Line weight here is a simply 1 00195 *ap_weight = 1.; 00196 // That's all 00197 return 0; 00198 } 00199 00200 // ===================================================================== 00201 // --------------------------------------------------------------------- 00202 // --------------------------------------------------------------------- 00203 // ===================================================================== 00204 00205 int iOptimizerOriginalAML::DataSpaceSpecificOperations( FLTNB a_data, FLTNB a_forwardModel, FLTNB* ap_backwardValues, 00206 FLTNB a_multiplicativeCorrections, FLTNB a_additiveCorrections, 00207 FLTNB a_quantificationFactor, oProjectionLine* ap_Line ) 00208 { 00209 // Implementation of the standard AML (bound is negative or null) 00210 if (m_bound<=0.) 00211 { 00212 // Compute the bound projection 00213 FLTNB shift = ForwardProject(ap_Line) * m_bound; 00214 // Compute denominator (this is the shifted model) 00215 FLTNB denominator = a_forwardModel - shift; 00216 // If the denominator is to close to zero, then ignore this data 00217 if (denominator<=m_dataSpaceDenominatorThreshold) *ap_backwardValues = 1.; 00218 // Otherwise, do it normally 00219 else 00220 { 00221 // Compute numerator (this is the shifted data) 00222 FLTNB numerator = a_data - shift; 00223 // Truncate to 0 if negative 00224 if (numerator<0.) *ap_backwardValues = 0.; 00225 // Otherwise, update backward values 00226 else *ap_backwardValues = numerator / denominator; 00227 } 00228 } 00229 // Implementation of AML with the bound taken to minus infinity (when the provided bound is positive) 00230 else 00231 { 00232 *ap_backwardValues = (a_data - a_forwardModel) / ForwardProject(ap_Line); 00233 } 00234 // That's all 00235 return 0; 00236 } 00237 00238 // ===================================================================== 00239 // --------------------------------------------------------------------- 00240 // --------------------------------------------------------------------- 00241 // ===================================================================== 00242 00243 int iOptimizerOriginalAML::ImageSpaceSpecificOperations( FLTNB a_currentImageValue, FLTNB* ap_newImageValue, 00244 FLTNB a_sensitivity, FLTNB* ap_correctionValues ) 00245 { 00246 // Compute image update factor 00247 FLTNB image_update_factor = *ap_correctionValues / a_sensitivity; 00248 // Update image for the standard AML (bound is negative or null) 00249 if (m_bound<=0.) *ap_newImageValue = m_bound + (a_currentImageValue - m_bound) * image_update_factor; 00250 // Update image for a bound taken to minus infinity (when the provided bound is positive) 00251 else *ap_newImageValue = a_currentImageValue + image_update_factor; 00252 // End 00253 return 0; 00254 } 00255 00256 // ===================================================================== 00257 // --------------------------------------------------------------------- 00258 // --------------------------------------------------------------------- 00259 // =====================================================================