![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class oComputeProjection 00004 00005 - separators: X 00006 - doxygen: X 00007 - default initialization: X 00008 - CASTOR_DEBUG: none 00009 - CASTOR_VERBOSE: X 00010 */ 00011 00019 #include "oComputeProjection.hh" 00020 00021 00022 // ===================================================================== 00023 // --------------------------------------------------------------------- 00024 // --------------------------------------------------------------------- 00025 // ===================================================================== 00026 /* 00027 \brief oComputeProjection constructor. 00028 \param ap_ImageDimensionsAndQuantification 00029 Initialize the member variables to their default values. 00030 */ 00031 oComputeProjection::oComputeProjection(oImageDimensionsAndQuantification* ap_ImageDimensionsAndQuantification) 00032 { 00033 mp_ID = ap_ImageDimensionsAndQuantification; 00034 m_noiseModelEnabled = false; 00035 } 00036 00037 00038 00039 00040 // ===================================================================== 00041 // --------------------------------------------------------------------- 00042 // --------------------------------------------------------------------- 00043 // ===================================================================== 00044 /* 00045 \brief oComputeProjection destructor. 00046 */ 00047 oComputeProjection::~oComputeProjection() 00048 { 00049 00050 } 00051 00052 00053 00054 00055 // ===================================================================== 00056 // --------------------------------------------------------------------- 00057 // --------------------------------------------------------------------- 00058 // ===================================================================== 00059 /* 00060 \fn DataUpdateStep 00061 \param ap_DataFile : pointer to the datafile which manage output data writing 00062 \param a2p_Line : pointer to the structure containing system matrix components 00063 \param ap_Image : pointer to the image space to access forward/attenuation/projection images 00064 \param ap_Event : pointer to the structure containing the event geometric indices 00065 \param a_fr : frame index 00066 \param a_rg : respiratory gate index 00067 \param a_cg : cardiac gate index 00068 \param a_timestamp : event timestamp 00069 \brief Perform the forward-projection, and send the result to the datafile for writing. 00070 Implements random Poisson noise. 00071 Write the result in the projection image matrix (only implemented for SPECT). 00072 \return 0 if success, positive value otherwise 00073 */ 00074 int oComputeProjection::DataUpdateStep(vDataFile* ap_DataFile, 00075 oProjectionLine* a2p_Line, 00076 oImageSpace* ap_Image, 00077 vEvent* ap_Event, 00078 int a_fr, 00079 int a_rg, 00080 int a_cg, 00081 int th, 00082 uint32_t a_timestamp) 00083 { 00084 00085 #ifdef CASTOR_VERBOSE 00086 if(m_verbose >=4) Cout("oComputeProjection::DataUpdateStep ... " << endl); 00087 #endif 00088 00089 FLTNB fp = 0., fp_atn = 0.; 00090 00091 // Projection with attenuation for PET 00092 //TODO : add fonction in oProjectionLine for PET attenuation ? 00093 if(ap_DataFile->GetDataType() == TYPE_PET) 00094 { 00095 // Forward project input image 00096 fp = a2p_Line->ForwardProject(ap_Image->m4p_forwardImage[a_fr][a_rg][a_cg]); 00097 00098 // Forward project attenuation if atn map has been provided 00099 if(ap_Image->m4p_attenuation[a_fr][a_rg][a_cg] != NULL) 00100 fp_atn = a2p_Line->ForwardProject(ap_Image->m4p_attenuation[a_fr][a_rg][a_cg]) * 0.1; // 0.1 -> convert in mm-1 00101 00102 // Apply attenuation factor 00103 fp *= std::exp(-fp_atn); 00104 00105 // Store Atn correction factor in Event if PET (not recorded in datafile if no atn map provided by the user) 00106 (dynamic_cast<iEventPET*>(ap_Event))->SetAttenuationCorrectionFactor(1/exp(-fp_atn)); 00107 } 00108 // Project with attenuation for SPECT 00109 else if (ap_DataFile->GetDataType() == TYPE_SPECT) 00110 fp = a2p_Line->ForwardProjectWithSPECTAttenuation(ap_Image->m4p_attenuation[a_fr][a_rg][a_cg] , 00111 ap_Image->m4p_forwardImage[a_fr][a_rg][a_cg] ); 00112 00113 00114 #ifdef CASTOR_VERBOSE 00115 if(m_verbose >=4) Cout("oComputeProjection::DataUpdateStep : Projeted value = "<< fp << endl); 00116 #endif 00117 00118 // Poisson noise 00119 if(m_noiseModelEnabled) 00120 { 00121 fp = GetPoissonNoise(fp); 00122 00123 #ifdef CASTOR_VERBOSE 00124 if(m_verbose >=4) Cout("oComputeProjection::DataUpdateStep : Projeted value with Poisson noise = "<< fp << endl); 00125 #endif 00126 } 00127 00128 // Write event in datafile 00129 ap_Event->SetEventValue(0, fp); 00130 ap_Event->SetTimeInMs(a_timestamp); 00131 00132 if(ap_DataFile->PROJ_WriteEvent(ap_Event, th) ) 00133 { 00134 Cerr("*****oComputeProjection::DataUpdateStep()-> Error while trying to write a projeted event" << endl); 00135 return 1; 00136 } 00137 00138 // Add the result to the projection image (only implemented for SPECT) 00139 if(ap_DataFile->GetDataType() == TYPE_SPECT) 00140 ap_Image->m2p_projectionImage[ap_Event->GetID1(0)][ap_Event->GetID2(0)] += fp; 00141 00142 return 0; 00143 } 00144 00145 00146 00147 00148 // ===================================================================== 00149 // --------------------------------------------------------------------- 00150 // --------------------------------------------------------------------- 00151 // ===================================================================== 00152 /* 00153 \fn InitNoiseModel 00154 \param aNoiseModel : string corresponding to a noise model 00155 \brief This is a premature implementation of noise model initialization 00156 for analytic simulator. 00157 Currently, only the Poisson noise can be selected. 00158 It will be replaced with a more efficient system when a more 00159 advanced analytical simulator will be implemented 00160 \return 0 if success, positive value otherwise 00161 */ 00162 int oComputeProjection::InitNoiseModel(string aNoiseModel) 00163 { 00164 if(m_verbose >=2) Cout("oComputeProjection::InitNoiseModel ... " << endl); 00165 00166 // Get the noise model name in the options and isolate the real model's options 00167 00168 // First check emptyness of the options 00169 if (aNoiseModel=="") 00170 { 00171 m_noiseModelEnabled = false; 00172 return 0; 00173 } 00174 else 00175 { 00176 string name_noise_model = ""; 00177 string list_options_noise_model = ""; 00178 string file_options_noise_model = ""; 00179 00180 // Search for a colon ":", this indicates that a configuration file is provided after the optimizer name 00181 size_t colon = aNoiseModel.find_first_of(":"); 00182 size_t comma = aNoiseModel.find_first_of(","); 00183 00184 // Case 1: we have a colon 00185 if (colon!=string::npos) 00186 { 00187 // Get the optimizer name before the colon 00188 name_noise_model = aNoiseModel.substr(0,colon); 00189 // Get the configuration file after the colon 00190 file_options_noise_model = aNoiseModel.substr(colon+1); 00191 // List of options is empty 00192 list_options_noise_model = ""; 00193 } 00194 // Case 2: we have a comma 00195 else if (comma!=string::npos) 00196 { 00197 // Get the optimizer name before the first comma 00198 name_noise_model = aNoiseModel.substr(0,comma); 00199 // Get the list of options after the first comma 00200 list_options_noise_model = aNoiseModel.substr(comma+1); 00201 // Configuration file is empty 00202 file_options_noise_model = ""; 00203 } 00204 // Case 3: no colon and no comma (a single model name) 00205 else 00206 { 00207 // Get the optimizer name 00208 name_noise_model = aNoiseModel; 00209 // Configuration file is empty 00210 file_options_noise_model = ""; 00211 // List of options is empty 00212 list_options_noise_model = ""; 00213 } 00214 00215 // Check if the noise model is known 00216 if (aNoiseModel == "poisson") 00217 { 00218 m_noiseModelEnabled = true; 00219 return 0; 00220 } 00221 else 00222 { 00223 Cerr("*****oComputeProjection::InitNoiseModel()-> Error while trying to initialize noise model. Model '"<<aNoiseModel<<"' is unknown"<<endl;); 00224 return 1; 00225 } 00226 } 00227 00228 return 0; 00229 } 00230 00231 00232 00233 00234 // ===================================================================== 00235 // --------------------------------------------------------------------- 00236 // --------------------------------------------------------------------- 00237 // ===================================================================== 00238 /* 00239 \fn GetPoissonNoise 00240 \param a_lambda : event rate 00241 \brief Generate a Poisson random variable from a provided lambda 00242 \details Depending on lambda, use either : 00243 - Poisson noise if a_lambda <= 60 00244 - Normal approximation otherwise 00245 \return a random value 00246 \todo Check Implementation 00247 \todo Use C++11 Poisson noise distribution ? 00248 */ 00249 uint32_t oComputeProjection::GetPoissonNoise(uint32_t a_lambda) 00250 { 00251 #ifdef CASTOR_VERBOSE 00252 if(m_verbose >=4) Cout("oComputeProjection::GetPoissonNoise ... " << endl); 00253 #endif 00254 00255 // Get instance of Random Number Generator 00256 sRNG* p_RNG = sRNG::GetInstance(); 00257 00258 double return_value = 0.; 00259 00260 // Normal approximation for projection with high value 00261 if (a_lambda > 60) 00262 { 00263 double rdm_1, rdm_2, w; 00264 00265 do 00266 { 00267 rdm_1 = 2.0 * p_RNG->GenerateRdmNber() - 1.0; 00268 rdm_2 = 2.0 * p_RNG->GenerateRdmNber() - 1.0; 00269 w = rdm_1*rdm_1 + rdm_2*rdm_2; 00270 } while (w > 1.0); 00271 00272 w = sqrt( (-2.0 * log(w) ) / w ); 00273 00274 return_value = a_lambda + sqrt(a_lambda)*rdm_1*w; 00275 00276 return_value = (return_value<0) ? 0 : round(return_value); 00277 } 00278 // Poisson noise 00279 else 00280 { 00281 double L = exp(-(a_lambda)); 00282 int n = 0; 00283 double p = 1; 00284 00285 do 00286 { 00287 n++; 00288 // Generate uniform random number u in [0,1] and let p ← p × u. 00289 p *= p_RNG->GenerateRdmNber(); 00290 } while (p > L); 00291 00292 return_value = n-1; 00293 } 00294 00295 return (uint32_t) return_value; 00296 }