 |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class vOptimizer 00004 00005 - separators: X 00006 - doxygen: 00007 - default initialization: X 00008 - CASTOR_DEBUG: 00009 - CASTOR_VERBOSE: 00010 */ 00011 00018 #include "vOptimizer.hh" 00019 00020 // ===================================================================== 00021 // --------------------------------------------------------------------- 00022 // --------------------------------------------------------------------- 00023 // ===================================================================== 00024 00025 vOptimizer::vOptimizer() 00026 { 00027 // Affect default values 00028 mp_ImageDimensionsAndQuantification = NULL; 00029 m_nbBackwardImages = 0; 00030 m_nbTOFBins = 0; 00031 m_initialValue = 0.; 00032 m_listmodeCompatibility = false; 00033 m_histogramCompatibility = false; 00034 m2p_forwardValues = NULL; 00035 m3p_backwardValues = NULL; 00036 m_dataMode = MODE_UNKNOWN; 00037 m_dataType = TYPE_UNKNOWN; 00038 m_optimizerFOMFlag = false; 00039 m_optimizerImageStatFlag = false; 00040 m4p_FOMLogLikelihood = NULL; 00041 m4p_FOMNbBins = NULL; 00042 m4p_FOMRMSE = NULL; 00043 m4p_FOMNbData = NULL; 00044 m_verbose = 0; 00045 m2p_attenuationImage = NULL; 00046 mp_imageStatNbVox = NULL; 00047 mp_imageStatMin = NULL; 00048 mp_imageStatMax = NULL; 00049 mp_imageStatMean = NULL; 00050 mp_imageStatVariance = NULL; 00051 mp_correctionStatMean = NULL; 00052 mp_correctionStatVariance = NULL; 00053 //m_penaltyEnergyFunctionDerivativesOrder = 0; 00054 } 00055 00056 // ===================================================================== 00057 // --------------------------------------------------------------------- 00058 // --------------------------------------------------------------------- 00059 // ===================================================================== 00060 00061 vOptimizer::~vOptimizer() 00062 { 00063 // Free forward and backward buffers 00064 if (m2p_forwardValues) 00065 { 00066 for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++) 00067 if (m2p_forwardValues[th]) free(m2p_forwardValues[th]); 00068 free(m2p_forwardValues); 00069 } 00070 if (m3p_backwardValues) 00071 { 00072 for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsMax(); th++) 00073 { 00074 if (m3p_backwardValues[th]) 00075 { 00076 for (int tof=0; tof<m_nbTOFBins; tof++) if (m3p_backwardValues[th][tof]) free(m3p_backwardValues[th][tof]); 00077 free(m3p_backwardValues[th]); 00078 } 00079 } 00080 free(m3p_backwardValues); 00081 } 00082 // Free pointers to attenuation images for SPECT 00083 if (m2p_attenuationImage) free(m2p_attenuationImage); 00084 // Delete all image update statistics tables 00085 if (m_optimizerImageStatFlag) 00086 { 00087 if (mp_imageStatNbVox) free(mp_imageStatNbVox); 00088 if (mp_imageStatMin) free(mp_imageStatMin); 00089 if (mp_imageStatMax) free(mp_imageStatMax); 00090 if (mp_imageStatMean) free(mp_imageStatMean); 00091 if (mp_imageStatVariance) free(mp_imageStatVariance); 00092 if (mp_correctionStatMean) free(mp_correctionStatMean); 00093 if (mp_correctionStatVariance) free(mp_correctionStatVariance); 00094 } 00095 // Delete FOM tables 00096 if (m_optimizerFOMFlag) 00097 { 00098 // Delete the log-likelihood table 00099 if (m4p_FOMLogLikelihood) 00100 { 00101 for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++) 00102 { 00103 if (m4p_FOMLogLikelihood[fr]) 00104 { 00105 for (int rg=0; rg<mp_ImageDimensionsAndQuantification->GetNbRespGates(); rg++) 00106 { 00107 if (m4p_FOMLogLikelihood[fr][rg]) 00108 { 00109 for (int cg=0; cg<mp_ImageDimensionsAndQuantification->GetNbCardGates(); cg++) 00110 if (m4p_FOMLogLikelihood[fr][rg][cg]) free(m4p_FOMLogLikelihood[fr][rg][cg]); 00111 free(m4p_FOMLogLikelihood[fr][rg]); 00112 } 00113 } 00114 free(m4p_FOMLogLikelihood[fr]); 00115 } 00116 } 00117 free(m4p_FOMLogLikelihood); 00118 } 00119 // Delete the RMSE table 00120 if (m4p_FOMRMSE) 00121 { 00122 for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++) 00123 { 00124 if (m4p_FOMRMSE[fr]) 00125 { 00126 for (int rg=0; rg<mp_ImageDimensionsAndQuantification->GetNbRespGates(); rg++) 00127 { 00128 if (m4p_FOMRMSE[fr][rg]) 00129 { 00130 for (int cg=0; cg<mp_ImageDimensionsAndQuantification->GetNbCardGates(); cg++) 00131 if (m4p_FOMRMSE[fr][rg][cg]) free(m4p_FOMRMSE[fr][rg][cg]); 00132 free(m4p_FOMRMSE[fr][rg]); 00133 } 00134 } 00135 free(m4p_FOMRMSE[fr]); 00136 } 00137 } 00138 free(m4p_FOMRMSE); 00139 } 00140 // Delete the nbBins table 00141 if (m4p_FOMNbBins) 00142 { 00143 for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++) 00144 { 00145 if (m4p_FOMNbBins[fr]) 00146 { 00147 for (int rg=0; rg<mp_ImageDimensionsAndQuantification->GetNbRespGates(); rg++) 00148 { 00149 if (m4p_FOMNbBins[fr][rg]) 00150 { 00151 for (int cg=0; cg<mp_ImageDimensionsAndQuantification->GetNbCardGates(); cg++) 00152 if (m4p_FOMNbBins[fr][rg][cg]) free(m4p_FOMNbBins[fr][rg][cg]); 00153 free(m4p_FOMNbBins[fr][rg]); 00154 } 00155 } 00156 free(m4p_FOMNbBins[fr]); 00157 } 00158 } 00159 free(m4p_FOMNbBins); 00160 } 00161 // Delete the nbData table 00162 if (m4p_FOMNbData) 00163 { 00164 for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++) 00165 { 00166 if (m4p_FOMNbData[fr]) 00167 { 00168 for (int rg=0; rg<mp_ImageDimensionsAndQuantification->GetNbRespGates(); rg++) 00169 { 00170 if (m4p_FOMNbData[fr][rg]) 00171 { 00172 for (int cg=0; cg<mp_ImageDimensionsAndQuantification->GetNbCardGates(); cg++) 00173 if (m4p_FOMNbData[fr][rg][cg]) free(m4p_FOMNbData[fr][rg][cg]); 00174 free(m4p_FOMNbData[fr][rg]); 00175 } 00176 } 00177 free(m4p_FOMNbData[fr]); 00178 } 00179 } 00180 free(m4p_FOMNbData); 00181 } 00182 } 00183 } 00184 00185 // ===================================================================== 00186 // --------------------------------------------------------------------- 00187 // --------------------------------------------------------------------- 00188 // ===================================================================== 00189 00190 void vOptimizer::ShowHelp() 00191 { 00192 // Call the specific help function from the children 00193 ShowHelpSpecific(); 00194 // Then, say if the child optimizer is compatible with histogram and/or list-mode data 00195 if (m_listmodeCompatibility && m_histogramCompatibility) cout << "This optimizer is compatible with both histogram and list-mode data." << endl; 00196 else if (m_listmodeCompatibility) cout << "This optimizer is only compatible with list-mode data." << endl; 00197 else if (m_histogramCompatibility) cout << "This optimizer is only compatible with histogram data." << endl; 00198 else cout << "This optimizer is a total non-sense, not compatible with histogram nor list-mode data !!!" << endl; 00199 } 00200 00201 // ===================================================================== 00202 // --------------------------------------------------------------------- 00203 // --------------------------------------------------------------------- 00204 // ===================================================================== 00205 00206 int vOptimizer::CheckParameters() 00207 { 00208 // Check image dimensions 00209 if (mp_ImageDimensionsAndQuantification==NULL) 00210 { 00211 Cerr("***** vOptimizer::CheckParameters() -> oImageDimensionsAndQuantification is null !" << endl); 00212 return 1; 00213 } 00214 // Check data mode 00215 if (m_dataMode==MODE_UNKNOWN || (m_dataMode!=MODE_LIST && m_dataMode!=MODE_HISTOGRAM)) 00216 { 00217 Cerr("***** vOptimizer::CheckParameters() -> No or meaningless data mode provided !" << endl); 00218 return 1; 00219 } 00220 // Check data type 00221 if (m_dataType==TYPE_UNKNOWN || (m_dataType!=TYPE_PET && m_dataType!=TYPE_SPECT && m_dataType!=TYPE_TRANSMISSION)) 00222 { 00223 Cerr("***** vOptimizer::CheckParameters() -> No or meaningless data type provided !" << endl); 00224 return 1; 00225 } 00226 // Check verbose level 00227 if (m_verbose<0) 00228 { 00229 Cerr("***** vOptimizer::CheckParameters() -> Verbose level is negative !" << endl); 00230 return 1; 00231 } 00232 // Listmode incompatibility 00233 if (m_dataMode==MODE_LIST && !m_listmodeCompatibility) 00234 { 00235 Cerr("***** vOptimizer::CheckParameters() -> The selected optimizer is not compatible with listmode data !" << endl); 00236 return 1; 00237 } 00238 // Histogram incompatibility 00239 if (m_dataMode==MODE_HISTOGRAM && !m_histogramCompatibility) 00240 { 00241 Cerr("***** vOptimizer::CheckParameters() -> The selected optimizer is not compatible with histogram data !" << endl); 00242 return 1; 00243 } 00244 // Send an error if FOM computation is required and verbose is null 00245 if (m_optimizerFOMFlag && m_verbose==0) 00246 { 00247 Cerr("***** vOptimizer::CheckParameters() -> Asked to compute FOMs in the data-space whereas the verbose is not set !" << endl); 00248 return 1; 00249 } 00250 // Send an error if image statistics computation is required and verbose is null 00251 if (m_optimizerImageStatFlag && m_verbose==0) 00252 { 00253 Cerr("***** vOptimizer::CheckParameters() -> Asked to compute image update statistics whereas the verbose is not set !" << endl); 00254 return 1; 00255 } 00256 // Send an error if FOM computation with list-mode data (it has no sense) 00257 if (m_optimizerFOMFlag && m_dataMode==MODE_LIST) 00258 { 00259 Cerr("***** vOptimizer::CheckParameters() -> Computing FOMs in the data-space with list-mode data is not possible !" << endl); 00260 return 1; 00261 } 00262 // Call the CheckSpecificParameters function of the child 00263 if (CheckSpecificParameters()) 00264 { 00265 Cerr("***** vOptimizer::CheckParameters() -> A problem occured while checking parameters specific to the optimizer module !" << endl); 00266 return 1; 00267 } 00268 // Normal end 00269 return 0; 00270 } 00271 00272 // ===================================================================== 00273 // --------------------------------------------------------------------- 00274 // --------------------------------------------------------------------- 00275 // ===================================================================== 00276 00277 int vOptimizer::Initialize() 00278 { 00279 // Allocate forward buffers (as many as threads and as many as TOF bins) 00280 m2p_forwardValues = (FLTNB**)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()*sizeof(FLTNB*)); 00281 for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++) 00282 m2p_forwardValues[th] = (FLTNB*)malloc(m_nbTOFBins*sizeof(FLTNB)); 00283 00284 // Allocate backward buffers (as many as threads then as many as TOF bins and as many as backward images) 00285 m3p_backwardValues = (FLTNB***)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsMax()*sizeof(FLTNB**)); 00286 for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsMax(); th++) 00287 { 00288 m3p_backwardValues[th] = (FLTNB**)malloc(m_nbTOFBins*sizeof(FLTNB*)); 00289 for (int tof=0; tof<m_nbTOFBins; tof++) 00290 m3p_backwardValues[th][tof] = (FLTNB*)malloc(m_nbBackwardImages*sizeof(FLTNB)); 00291 } 00292 00293 // Allocate pointers to the current attenuation image for SPECT attenuation correction 00294 m2p_attenuationImage = (FLTNB**)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()*sizeof(FLTNB*)); 00295 for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++) m2p_attenuationImage[th] = NULL; 00296 00297 // Allocate the thread-safe tables for image update statistics computation 00298 if (m_optimizerImageStatFlag) 00299 { 00300 mp_imageStatNbVox = (INTNB*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation()*sizeof(INTNB)); 00301 mp_imageStatMin = (FLTNB*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation()*sizeof(FLTNB)); 00302 mp_imageStatMax = (FLTNB*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation()*sizeof(FLTNB)); 00303 mp_imageStatMean = (double*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation()*sizeof(double)); 00304 mp_imageStatVariance = (double*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation()*sizeof(double)); 00305 mp_correctionStatMean = (double*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation()*sizeof(double)); 00306 mp_correctionStatVariance = (double*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation()*sizeof(double)); 00307 } 00308 00309 // We allocate the tables that hold the figures-of-merit if FOM flag is on 00310 if (m_optimizerFOMFlag) 00311 { 00312 m4p_FOMNbBins = (uint64_t****)malloc(mp_ImageDimensionsAndQuantification->GetNbTimeFrames()*sizeof(uint64_t***)); 00313 m4p_FOMNbData = (double****)malloc(mp_ImageDimensionsAndQuantification->GetNbTimeFrames()*sizeof(double***)); 00314 m4p_FOMLogLikelihood = (FLTNB****)malloc(mp_ImageDimensionsAndQuantification->GetNbTimeFrames()*sizeof(FLTNB***)); 00315 m4p_FOMRMSE = (FLTNB****)malloc(mp_ImageDimensionsAndQuantification->GetNbTimeFrames()*sizeof(FLTNB***)); 00316 for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++) 00317 { 00318 m4p_FOMNbBins[fr] = (uint64_t***)malloc(mp_ImageDimensionsAndQuantification->GetNbRespGates()*sizeof(uint64_t**)); 00319 m4p_FOMNbData[fr] = (double***)malloc(mp_ImageDimensionsAndQuantification->GetNbRespGates()*sizeof(double**)); 00320 m4p_FOMLogLikelihood[fr] = (FLTNB***)malloc(mp_ImageDimensionsAndQuantification->GetNbRespGates()*sizeof(FLTNB**)); 00321 m4p_FOMRMSE[fr] = (FLTNB***)malloc(mp_ImageDimensionsAndQuantification->GetNbRespGates()*sizeof(FLTNB**)); 00322 for (int rg=0; rg<mp_ImageDimensionsAndQuantification->GetNbRespGates(); rg++) 00323 { 00324 m4p_FOMNbBins[fr][rg] = (uint64_t**)malloc(mp_ImageDimensionsAndQuantification->GetNbCardGates()*sizeof(uint64_t*)); 00325 m4p_FOMNbData[fr][rg] = (double**)malloc(mp_ImageDimensionsAndQuantification->GetNbCardGates()*sizeof(double*)); 00326 m4p_FOMLogLikelihood[fr][rg] = (FLTNB**)malloc(mp_ImageDimensionsAndQuantification->GetNbCardGates()*sizeof(FLTNB*)); 00327 m4p_FOMRMSE[fr][rg] = (FLTNB**)malloc(mp_ImageDimensionsAndQuantification->GetNbCardGates()*sizeof(FLTNB*)); 00328 for (int cg=0; cg<mp_ImageDimensionsAndQuantification->GetNbCardGates(); cg++) 00329 { 00330 m4p_FOMNbBins[fr][rg][cg] = (uint64_t*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()*sizeof(uint64_t)); 00331 m4p_FOMNbData[fr][rg][cg] = (double*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()*sizeof(double)); 00332 m4p_FOMLogLikelihood[fr][rg][cg] = (FLTNB*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()*sizeof(FLTNB)); 00333 m4p_FOMRMSE[fr][rg][cg] = (FLTNB*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()*sizeof(FLTNB)); 00334 } 00335 } 00336 } 00337 } 00338 00339 // Call the specific initialization function of the child 00340 if (InitializeSpecific()) 00341 { 00342 Cerr("***** vOptimizer::Initialize() -> A problem occured while initializing stuff specific to the optimizer module !" << endl); 00343 return 1; 00344 } 00345 00346 // Normal end 00347 return 0; 00348 } 00349 00350 // ===================================================================== 00351 // --------------------------------------------------------------------- 00352 // --------------------------------------------------------------------- 00353 // ===================================================================== 00354 00355 int vOptimizer::PreDataUpdateStep(int a_iteration, int a_nbIterations, int a_subset, int a_nbSubsets) 00356 { 00357 // Simply reset FOMs if activated 00358 if (m_optimizerFOMFlag) 00359 { 00360 for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++) 00361 for (int rg=0; rg<mp_ImageDimensionsAndQuantification->GetNbRespGates(); rg++) 00362 for (int cg=0; cg<mp_ImageDimensionsAndQuantification->GetNbCardGates(); cg++) 00363 for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++) 00364 { 00365 m4p_FOMNbBins[fr][rg][cg][th] = 0; 00366 m4p_FOMNbData[fr][rg][cg][th] = 0.; 00367 m4p_FOMLogLikelihood[fr][rg][cg][th] = 0.; 00368 m4p_FOMRMSE[fr][rg][cg][th] = 0.; 00369 } 00370 } 00371 // Then call the specific pre-update-step function from the child optimizer 00372 if (PreDataUpdateSpecificStep(a_iteration, a_nbIterations, a_subset, a_nbSubsets)) 00373 { 00374 Cerr("***** vOptimizer::PreDataUpdateStep() -> A problem occured while applying the specific pre-data-update step of the optimizer module !" << endl); 00375 return 1; 00376 } 00377 // Normal end 00378 return 0; 00379 } 00380 00381 // ===================================================================== 00382 // --------------------------------------------------------------------- 00383 // --------------------------------------------------------------------- 00384 // ===================================================================== 00385 00386 int vOptimizer::PreDataUpdateSpecificStep(int a_iteration, int a_nbIterations, int a_subset, int a_nbSubsets) 00387 { 00388 // By default, just do nothing here. This function is designed to be overloaded by the specific optimizer if needed 00389 return 0; 00390 } 00391 00392 // ===================================================================== 00393 // --------------------------------------------------------------------- 00394 // --------------------------------------------------------------------- 00395 // ===================================================================== 00396 00397 int vOptimizer::PostDataUpdateStep(int a_iteration, int a_nbIterations, int a_subset, int a_nbSubsets) 00398 { 00399 // Print out the FOMs if activated 00400 if (m_optimizerFOMFlag) 00401 { 00402 // Verbose 00403 Cout("vOptimizer::PostDataUpdateStep() -> Data-space figures-of-merit" << endl); 00404 // Number of spaces for printing 00405 string spaces_fr = ""; if (mp_ImageDimensionsAndQuantification->GetNbTimeFrames()>1) spaces_fr += " "; 00406 string spaces_rg = ""; if (mp_ImageDimensionsAndQuantification->GetNbRespGates()>1) spaces_rg += " "; 00407 string spaces_cg = ""; if (mp_ImageDimensionsAndQuantification->GetNbCardGates()>1) spaces_cg += " "; 00408 // Loop on frames and gates 00409 for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++) 00410 { 00411 // Verbose 00412 if (mp_ImageDimensionsAndQuantification->GetNbTimeFrames()>1) 00413 Cout(spaces_fr << "Frame " << fr+1 << " on " << mp_ImageDimensionsAndQuantification->GetNbTimeFrames() << endl); 00414 for (int rg=0; rg<mp_ImageDimensionsAndQuantification->GetNbRespGates(); rg++) 00415 { 00416 // Verbose 00417 if (mp_ImageDimensionsAndQuantification->GetNbRespGates()>1) 00418 Cout(spaces_fr << spaces_rg << "Respiratory gate " << rg+1 << " on " << mp_ImageDimensionsAndQuantification->GetNbRespGates() << endl); 00419 for (int cg=0; cg<mp_ImageDimensionsAndQuantification->GetNbCardGates(); cg++) 00420 { 00421 // Verbose 00422 if (mp_ImageDimensionsAndQuantification->GetNbCardGates()>1) 00423 Cout(spaces_fr << spaces_rg << spaces_cg << "Cardiac gate " << cg+1 << " on " << mp_ImageDimensionsAndQuantification->GetNbCardGates() << endl); 00424 // Merge thread results 00425 for (int th=1; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++) 00426 { 00427 m4p_FOMNbBins[fr][rg][cg][0] += m4p_FOMNbBins[fr][rg][cg][th]; 00428 m4p_FOMNbData[fr][rg][cg][0] += m4p_FOMNbData[fr][rg][cg][th]; 00429 m4p_FOMLogLikelihood[fr][rg][cg][0] += m4p_FOMLogLikelihood[fr][rg][cg][th]; 00430 m4p_FOMRMSE[fr][rg][cg][0] += m4p_FOMRMSE[fr][rg][cg][th]; 00431 } 00432 // Finish mean number of data counts per channel 00433 m4p_FOMNbData[fr][rg][cg][0] /= ((double)(m4p_FOMNbBins[fr][rg][cg][0])); 00434 // Finish RMSE computation 00435 m4p_FOMRMSE[fr][rg][cg][0] /= ((FLTNB)(m4p_FOMNbBins[fr][rg][cg][0])); 00436 m4p_FOMRMSE[fr][rg][cg][0] = sqrt(m4p_FOMRMSE[fr][rg][cg][0]); 00437 // Print out 00438 Cout(spaces_fr << spaces_rg << spaces_cg << " --> Number of data channels used in optimization: " << m4p_FOMNbBins[fr][rg][cg][0] << endl); 00439 Cout(spaces_fr << spaces_rg << spaces_cg << " --> Mean number of data counts per channel: " << m4p_FOMNbData[fr][rg][cg][0] << endl); 00440 Cout(spaces_fr << spaces_rg << spaces_cg << " --> Log-likelihood: " << m4p_FOMLogLikelihood[fr][rg][cg][0] << endl); 00441 Cout(spaces_fr << spaces_rg << spaces_cg << " --> RMSE: " << m4p_FOMRMSE[fr][rg][cg][0] << endl); 00442 } 00443 } 00444 } 00445 } 00446 // Then call the specific post-update-step function from the child optimizer 00447 if (PostDataUpdateSpecificStep(a_iteration, a_nbIterations, a_subset, a_nbSubsets)) 00448 { 00449 Cerr("***** vOptimizer::PostDataUpdateStep() -> A problem occured while applying the specific post-data-update step of the optimizer module !" << endl); 00450 return 1; 00451 } 00452 // Normal end 00453 return 0; 00454 } 00455 00456 // ===================================================================== 00457 // --------------------------------------------------------------------- 00458 // --------------------------------------------------------------------- 00459 // ===================================================================== 00460 00461 int vOptimizer::PostDataUpdateSpecificStep(int a_iteration, int a_nbIterations, int a_subset, int a_nbSubsets) 00462 { 00463 // By default, just do nothing here. This function is designed to be overloaded by the specific optimizer if needed 00464 return 0; 00465 } 00466 00467 // ===================================================================== 00468 // --------------------------------------------------------------------- 00469 // --------------------------------------------------------------------- 00470 // ===================================================================== 00471 00472 int vOptimizer::DataStep1ForwardProjectModel( oProjectionLine* ap_Line, oImageSpace* ap_Image, vEvent* ap_Event, 00473 int a_bed, int a_timeFrame, int a_respGate, int a_cardGate, 00474 int a_th ) 00475 { 00476 // ======================================================================= 00477 // 4D forward projection including framing, respiratory and cardiac gating 00478 // ======================================================================= 00479 00480 // Clean buffers 00481 for (int b=0; b<ap_Line->GetNbTOFBins(); b++) 00482 { 00483 m2p_forwardValues[a_th][b] = 0.; 00484 for (int img=0; img<m_nbBackwardImages; img++) m3p_backwardValues[a_th][b][img] = 0.; 00485 } 00486 00487 // Loop over time basis functions 00488 for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++) 00489 { 00490 // Get frame/basis coefficient and continue if null 00491 FLTNB time_basis_coef = mp_ImageDimensionsAndQuantification->GetTimeBasisCoefficient(tbf,a_timeFrame); 00492 if (time_basis_coef==0.) continue; 00493 // Loop over respiratory basis functions 00494 for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++) 00495 { 00496 // Get resp_gate/basis coefficient and continue if null 00497 FLTNB resp_basis_coef = mp_ImageDimensionsAndQuantification->GetRespBasisCoefficient(rbf,a_respGate); 00498 if (resp_basis_coef==0.) continue; 00499 // Loop over cardiac basis functions 00500 for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++) 00501 { 00502 // Get card_gate_basis coefficient and continue if null 00503 FLTNB card_basis_coef = mp_ImageDimensionsAndQuantification->GetCardBasisCoefficient(cbf,a_cardGate); 00504 if (card_basis_coef==0.) continue; 00505 // Compute global coefficient 00506 FLTNB global_basis_coef = time_basis_coef * resp_basis_coef * card_basis_coef; 00507 // Loop over TOF bins 00508 for (int b=0; b<ap_Line->GetNbTOFBins(); b++) 00509 { 00510 // Set the current TOF bin 00511 ap_Line->SetCurrentTOFBin(b); 00512 // Project the image and apply dynamic coefficient converions 00513 m2p_forwardValues[a_th][b] += ForwardProject(ap_Line,ap_Image->m4p_forwardImage[tbf][rbf][cbf]) * global_basis_coef; 00514 } 00515 } 00516 } 00517 } 00518 00519 // ===================================== 00520 // Include the additive correction terms 00521 // ===================================== 00522 00523 // Add additive terms (we multiply by the frame duration because the additive terms are provided as rates) 00524 for (int b=0; b<ap_Line->GetNbTOFBins(); b++) m2p_forwardValues[a_th][b] += ap_Event->GetAdditiveCorrections(b) * mp_ImageDimensionsAndQuantification->GetFrameDurationInSec(a_bed, a_timeFrame); 00525 00526 // End 00527 return 0; 00528 } 00529 00530 // ===================================================================== 00531 // --------------------------------------------------------------------- 00532 // --------------------------------------------------------------------- 00533 // ===================================================================== 00534 00535 int vOptimizer::DataStep2Optional( oProjectionLine* ap_Line, oImageSpace* ap_Image, vEvent* ap_Event, 00536 int a_bed, int a_timeFrame, int a_respGate, int a_cardGate, 00537 int a_iteration, int a_th ) 00538 { 00539 // Deliberately do nothing! 00540 return 0; 00541 } 00542 00543 // ===================================================================== 00544 // --------------------------------------------------------------------- 00545 // --------------------------------------------------------------------- 00546 // ===================================================================== 00547 00548 int vOptimizer::DataStep3BackwardProjectSensitivity( oProjectionLine* ap_Line, oImageSpace* ap_Image, vEvent* ap_Event, 00549 int a_bed, int a_timeFrame, int a_respGate, int a_cardGate, 00550 int a_th ) 00551 { 00552 // Loop over TOF bins 00553 for (int b=0; b<ap_Line->GetNbTOFBins(); b++) 00554 { 00555 // Set the current TOF bin of the projection-line 00556 ap_Line->SetCurrentTOFBin(b); 00557 // The weight associated to the line to be backprojected 00558 FLTNB weight = 0.; 00559 // Call the function to compute specific sensitivity 00560 if (SensitivitySpecificOperations 00561 ( 00562 ap_Event->GetEventValue(b), // 1: the measured data 00563 m2p_forwardValues[a_th][b], // 2: the forward model 00564 &weight, // 3: the sensitivity weight 00565 ap_Event->GetMultiplicativeCorrections(), // 4: the multiplicative corrections 00566 ap_Event->GetAdditiveCorrections(b)*mp_ImageDimensionsAndQuantification->GetFrameDurationInSec(a_bed, a_timeFrame), // 5: the additive corrections 00567 mp_ImageDimensionsAndQuantification->GetQuantificationFactor(a_bed,a_timeFrame, a_respGate, a_cardGate), // 6: the quantification factor 00568 ap_Line // 7: the projection line 00569 )) 00570 { 00571 Cerr("***** vOptimizer::DataStep3BackwardProjectSensitivity() -> A problem occured while performing the sensitivity specific operations !" << endl); 00572 return 1; 00573 } 00574 // Back-project the weight into the sensitivity matrix 00575 BackwardProject(ap_Line, ap_Image->m5p_sensitivity[a_th][a_timeFrame][a_respGate][a_cardGate], weight); 00576 } 00577 00578 // End 00579 return 0; 00580 } 00581 00582 // ===================================================================== 00583 // --------------------------------------------------------------------- 00584 // --------------------------------------------------------------------- 00585 // ===================================================================== 00586 00587 int vOptimizer::DataStep4Optional( oProjectionLine* ap_Line, oImageSpace* ap_Image, vEvent* ap_Event, 00588 int a_bed, int a_timeFrame, int a_respGate, int a_cardGate, 00589 int a_iteration, int a_th ) 00590 { 00591 // Deliberately do nothing! 00592 return 0; 00593 } 00594 00595 // ===================================================================== 00596 // --------------------------------------------------------------------- 00597 // --------------------------------------------------------------------- 00598 // ===================================================================== 00599 00600 int vOptimizer::DataStep5ComputeCorrections( oProjectionLine* ap_Line, vEvent* ap_Event, 00601 int a_bed, int a_timeFrame, int a_respGate, int a_cardGate, 00602 int a_th ) 00603 { 00604 // Loop on TOF bins 00605 for (int b=0; b<ap_Line->GetNbTOFBins(); b++) 00606 { 00607 // Set the current TOF bin of the projection-line 00608 ap_Line->SetCurrentTOFBin(b); 00609 // Then call the pure virtual function for specific data space operations 00610 if (DataSpaceSpecificOperations 00611 ( 00612 ap_Event->GetEventValue(b), // 1: the measured data 00613 m2p_forwardValues[a_th][b], // 2: the forward model 00614 m3p_backwardValues[a_th][b], // 3: the backward values (the result) 00615 ap_Event->GetMultiplicativeCorrections(), // 4: the multiplicative corrections 00616 ap_Event->GetAdditiveCorrections(b)*mp_ImageDimensionsAndQuantification->GetFrameDurationInSec(a_bed, a_timeFrame), // 5: the additive corrections 00617 mp_ImageDimensionsAndQuantification->GetQuantificationFactor(a_bed,a_timeFrame, a_respGate, a_cardGate), // 6: the quantification factor 00618 ap_Line // 7: the projection line 00619 )) 00620 { 00621 Cerr("***** vOptimizer::DataStep5ComputeCorrections() -> A problem occured while performing specific data space operations !" << endl); 00622 return 1; 00623 } 00624 } 00625 00626 // End 00627 return 0; 00628 } 00629 00630 // ===================================================================== 00631 // --------------------------------------------------------------------- 00632 // --------------------------------------------------------------------- 00633 // ===================================================================== 00634 00635 int vOptimizer::DataStep6Optional( oProjectionLine* ap_Line, oImageSpace* ap_Image, vEvent* ap_Event, 00636 int a_bed, int a_timeFrame, int a_respGate, int a_cardGate, 00637 int a_iteration, int a_th ) 00638 { 00639 // Deliberately do nothing! 00640 return 0; 00641 } 00642 00643 // ===================================================================== 00644 // --------------------------------------------------------------------- 00645 // --------------------------------------------------------------------- 00646 // ===================================================================== 00647 00648 int vOptimizer::DataStep7BackwardProjectCorrections( oProjectionLine* ap_Line, oImageSpace* ap_Image, vEvent* ap_Event, 00649 int a_bed, int a_timeFrame, int a_respGate, int a_cardGate, 00650 int a_th ) 00651 { 00652 // ======================================================================== 00653 // 4D backward projection including framing, respiratory and cardiac gating 00654 // ======================================================================== 00655 00656 // Loop over time basis functions 00657 for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++) 00658 { 00659 // Get frame/basis coefficient and continue if null 00660 FLTNB time_basis_coef = mp_ImageDimensionsAndQuantification->GetTimeBasisCoefficient(tbf,a_timeFrame); 00661 if (time_basis_coef==0.) continue; 00662 // Loop over respiratory basis functions 00663 for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++) 00664 { 00665 // Get resp_gate/basis coefficient and continue if null 00666 FLTNB resp_basis_coef = mp_ImageDimensionsAndQuantification->GetRespBasisCoefficient(rbf,a_respGate); 00667 if (resp_basis_coef==0.) continue; 00668 // Loop over cardiac basis functions 00669 for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++) 00670 { 00671 // Get card_gate_basis coefficient and continue if null 00672 FLTNB card_basis_coef = mp_ImageDimensionsAndQuantification->GetCardBasisCoefficient(cbf,a_cardGate); 00673 if (card_basis_coef==0.) continue; 00674 // Compute global coefficient 00675 FLTNB global_basis_coef = time_basis_coef * resp_basis_coef * card_basis_coef; 00676 // Loop over TOF bins 00677 for (int b=0; b<ap_Line->GetNbTOFBins(); b++) 00678 { 00679 // Set the current TOF bin of the projection-line 00680 ap_Line->SetCurrentTOFBin(b); 00681 // Loop over backprojection images 00682 for (int img=0; img<m_nbBackwardImages; img++) 00683 { 00684 // Continue if backproj is null or infinity 00685 if (m3p_backwardValues[a_th][b][img]==0. || !isfinite(m3p_backwardValues[a_th][b][img])) continue; 00686 // Backward project into the backward image 00687 BackwardProject( ap_Line, ap_Image->m6p_backwardImage[img][a_th][tbf][rbf][cbf] , m3p_backwardValues[a_th][b][img] * global_basis_coef ); 00688 } 00689 } 00690 } 00691 } 00692 } 00693 00694 // End 00695 return 0; 00696 } 00697 00698 // ===================================================================== 00699 // --------------------------------------------------------------------- 00700 // --------------------------------------------------------------------- 00701 // ===================================================================== 00702 00703 int vOptimizer::DataStep8ComputeFOM( oProjectionLine* ap_Line, vEvent* ap_Event, 00704 int a_timeFrame, int a_respGate, int a_cardGate, 00705 int a_thread ) 00706 { 00707 // Loop on TOF bins 00708 for (int b=0; b<ap_Line->GetNbTOFBins(); b++) 00709 { 00710 // Get the model and the data 00711 FLTNB model = m2p_forwardValues[a_thread][b]; 00712 FLTNB data = ap_Event->GetEventValue(b); 00713 // Update log likelihood only if model is strictly positive (to avoid numerical issues, we use a threshold at e-10) 00714 if (model>1.e-10) m4p_FOMLogLikelihood[a_timeFrame][a_respGate][a_cardGate][a_thread] += data*log(model)-model; 00715 // Update RMSE 00716 m4p_FOMRMSE[a_timeFrame][a_respGate][a_cardGate][a_thread] += (data-model)*(data-model); 00717 // Update number of data channels 00718 m4p_FOMNbBins[a_timeFrame][a_respGate][a_cardGate][a_thread]++; 00719 // Update the number of data counts 00720 m4p_FOMNbData[a_timeFrame][a_respGate][a_cardGate][a_thread] += ((double)data); 00721 } 00722 // Normal end 00723 return 0; 00724 } 00725 00726 // ===================================================================== 00727 // --------------------------------------------------------------------- 00728 // --------------------------------------------------------------------- 00729 // ===================================================================== 00730 00731 int vOptimizer::UpdateVisitedVoxels( oImageSpace* ap_Image ) 00732 { 00733 // Verbose 00734 if (m_verbose>=3) Cout("vOptimizer::UpdateVisitedVoxels() -> Tick visited voxels based on sensitivity" << endl); 00735 00736 // We loop on the sensitivity for the first frame/gates and look if voxels were 'visited' by LORs 00737 int thread0 = 0; int frame0 = 0; int resp_gate0 = 0; int card_gate0 = 0; 00738 for (int v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++) 00739 // If the voxel has been visited, we add 1. to the ap_Image->mp_visitedVoxelsImage 00740 if (ap_Image->m5p_sensitivity[thread0][frame0][resp_gate0][card_gate0][v]!=0.) ap_Image->mp_visitedVoxelsImage[v] += 1.; 00741 00742 // End 00743 return 0; 00744 } 00745 00746 // ===================================================================== 00747 // --------------------------------------------------------------------- 00748 // --------------------------------------------------------------------- 00749 // ===================================================================== 00750 00751 int vOptimizer::ImageUpdateStep( oImageSpace* ap_Image, int a_iteration, int a_nbSubsets ) 00752 { 00753 // Verbose 00754 if (m_verbose>=2 || m_optimizerImageStatFlag) Cout("vOptimizer::ImageUpdateStep() -> Start image update" << endl); 00755 // Number of spaces for nice printing 00756 string spaces_tbf = ""; if (mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions()>1) spaces_tbf += " "; 00757 string spaces_rbf = ""; if (mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions()>1) spaces_rbf += " "; 00758 string spaces_cbf = ""; if (mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions()>1) spaces_cbf += " "; 00759 // Loop over time basis functions 00760 for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++) 00761 { 00762 // Verbose 00763 if (mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions()>1 && (m_verbose>=3 || m_optimizerImageStatFlag)) 00764 { 00765 if (mp_ImageDimensionsAndQuantification->GetTimeStaticFlag()) Cout(spaces_tbf << "--> Time frame " << tbf+1 << endl); 00766 else Cout(spaces_tbf << "--> Time basis function: " << tbf+1 << endl); 00767 } 00768 // Loop over respiratory basis functions 00769 for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++) 00770 { 00771 // Verbose 00772 if (mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions()>1 && (m_verbose>=3 || m_optimizerImageStatFlag)) 00773 { 00774 if (mp_ImageDimensionsAndQuantification->GetRespStaticFlag()) Cout(spaces_tbf << spaces_rbf << "--> Respiratory gate " << rbf+1 << endl); 00775 else Cout(spaces_tbf << spaces_rbf << "--> Respiratory basis function: " << rbf+1 << endl); 00776 } 00777 // Loop over cardiac basis functions 00778 for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++) 00779 { 00780 // Verbose 00781 if (mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions()>1 && (m_verbose>=3 || m_optimizerImageStatFlag)) 00782 { 00783 if (mp_ImageDimensionsAndQuantification->GetCardStaticFlag()) Cout(spaces_tbf << spaces_rbf << spaces_cbf << "--> Cardiac gate " << cbf+1 << endl); 00784 else Cout(spaces_tbf << spaces_rbf << spaces_cbf << "--> Cardiac basis function: " << cbf+1 << endl); 00785 } 00786 // Set the number of threads for image computation 00787 #ifdef CASTOR_OMP 00788 omp_set_num_threads(mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation()); 00789 #endif 00790 // Reset counter for image stats 00791 if (m_optimizerImageStatFlag) 00792 { 00793 for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation(); th++) 00794 { 00795 mp_imageStatNbVox[th] = 0; 00796 mp_imageStatMin[th] = ap_Image->m4p_image[tbf][rbf][cbf][0]; 00797 mp_imageStatMax[th] = ap_Image->m4p_image[tbf][rbf][cbf][0]; 00798 mp_imageStatMean[th] = 0.; 00799 mp_imageStatVariance[th] = 0.; 00800 mp_correctionStatMean[th] = 0.; 00801 mp_correctionStatVariance[th] = 0.; 00802 } 00803 } 00804 // OMP loop over voxels 00805 int v; 00806 bool error = false; 00807 #pragma omp parallel for private(v) schedule(guided) 00808 for (v=0 ; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ() ; v++) 00809 { 00810 // Get the thread index 00811 int th = 0; 00812 #ifdef CASTOR_OMP 00813 th = omp_get_thread_num(); 00814 #endif 00815 // Compute sensitivity 00816 FLTNB sensitivity = ComputeSensitivity(ap_Image->m5p_sensitivity[0],tbf,rbf,cbf,v,a_nbSubsets); 00817 // If the sensitivity is negative or null, we skip this voxel 00818 if (sensitivity<=0.) continue; 00819 // Keep the current image value in a buffer for update statistics 00820 FLTNB old_image_value = ap_Image->m4p_image[tbf][rbf][cbf][v]; 00821 // Fill class-local buffer of backward values (thread-safe) 00822 for (int nb=0; nb<m_nbBackwardImages; nb++) 00823 m3p_backwardValues[th][0][nb] = ap_Image->m6p_backwardImage[nb][0][tbf][rbf][cbf][v]; 00824 // Then call the pure virtual function for specific data space operations 00825 if (ImageSpaceSpecificOperations 00826 ( 00827 ap_Image->m4p_image[tbf][rbf][cbf][v], 00828 // Note: have to be carefull here because we pass rbf and cbf as the actual respgate and cardgate. It may be wrong if one would like 00829 // to perform MLAA reconstruction but using PET basis functions different from diracs for example... In this case, we would get 00830 // segmentation faults. But actually, this function would be overloaded by MLAA anyway... 00831 //ap_Image->m4p_attenuation != NULL ? ap_Image->m4p_attenuation[tbf][rbf][cbf][v] : 0. , 00832 &ap_Image->m4p_image[tbf][rbf][cbf][v], 00833 //ap_Image->m4p_attenuation != NULL ? &ap_Image->m4p_attenuation[tbf][rbf][cbf][v] : NULL , 00834 sensitivity, 00835 m3p_backwardValues[th][0] 00836 )) 00837 { 00838 Cerr("***** vOptimizer::ImageUpdateStep() -> A problem occured in while performing image space specific operations of the optimizer !" << endl); 00839 error = true; 00840 } 00841 // Do some image statistics 00842 if (m_optimizerImageStatFlag) 00843 { 00844 if (mp_imageStatMin[th]>ap_Image->m4p_image[tbf][rbf][cbf][v]) mp_imageStatMin[th] = ap_Image->m4p_image[tbf][rbf][cbf][v]; 00845 if (mp_imageStatMax[th]<ap_Image->m4p_image[tbf][rbf][cbf][v]) mp_imageStatMax[th] = ap_Image->m4p_image[tbf][rbf][cbf][v]; 00846 // The one pass algorithm is used to compute the variance here for both the image and correction step. 00847 // Do not change anyhting to the order of the operations ! 00848 mp_imageStatNbVox[th]++; 00849 double sample = ((double)(ap_Image->m4p_image[tbf][rbf][cbf][v])); 00850 double delta = sample - mp_imageStatMean[th]; 00851 mp_imageStatMean[th] += delta / ((double)(mp_imageStatNbVox[th])); 00852 mp_imageStatVariance[th] += delta * (sample - mp_imageStatMean[th]); 00853 // Same one-pass variance algorithm for the correction step 00854 sample = ((double)(ap_Image->m4p_image[tbf][rbf][cbf][v] - old_image_value)); 00855 delta = sample - mp_correctionStatMean[th]; 00856 mp_correctionStatMean[th] += delta / ((double)(mp_imageStatNbVox[th])); 00857 mp_correctionStatVariance[th] += delta * (sample - mp_correctionStatMean[th]); 00858 } 00859 } 00860 // Check for errors (we cannot stop the multi-threaded loop so we do it here) 00861 if (error) return 1; 00862 // Finish image update statistics computations: these operations are specific to the merging of the 00863 // multi-threaded one-pass variance estimations used above. 00864 if (m_optimizerImageStatFlag) 00865 { 00866 for (int th=1; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation(); th++) 00867 { 00868 // Image minimum and maximum 00869 if (mp_imageStatMin[0]>mp_imageStatMin[th]) mp_imageStatMin[0] = mp_imageStatMin[th]; 00870 if (mp_imageStatMax[0]<mp_imageStatMax[th]) mp_imageStatMax[0] = mp_imageStatMax[th]; 00871 // Cast the counts into double 00872 double count1 = ((double)(mp_imageStatNbVox[0])); 00873 double count2 = ((double)(mp_imageStatNbVox[th])); 00874 // Update the count 00875 mp_imageStatNbVox[0] += mp_imageStatNbVox[th]; 00876 double count12 = ((double)(mp_imageStatNbVox[0])); 00877 // Compute the delta mean 00878 double delta = mp_imageStatMean[0] - mp_imageStatMean[th]; 00879 // Update the variance 00880 mp_imageStatVariance[0] += mp_imageStatVariance[th] + delta * delta * count1 * count2 / count12; 00881 // Update the mean 00882 mp_imageStatMean[0] = (count1*mp_imageStatMean[0] + count2*mp_imageStatMean[th]) / count12; 00883 // Do the same for the correction step 00884 delta = mp_correctionStatMean[0] - mp_correctionStatMean[th]; 00885 mp_correctionStatVariance[0] += mp_correctionStatVariance[th] + delta * delta * count1 * count2 / count12; 00886 mp_correctionStatMean[0] = (count1*mp_correctionStatMean[0] + count2*mp_correctionStatMean[th]) / count12; 00887 } 00888 // Finish variance computation 00889 mp_imageStatVariance[0] /= ((double)(mp_imageStatNbVox[0])); 00890 mp_correctionStatVariance[0] /= ((double)(mp_imageStatNbVox[0])); 00891 // Print out 00892 Cout(spaces_tbf << spaces_rbf << spaces_cbf << " --> Image update step " 00893 << " | mean: " << mp_correctionStatMean[0] 00894 << " | stdv: " << sqrt(mp_correctionStatVariance[0]) << endl); 00895 Cout(spaces_tbf << spaces_rbf << spaces_cbf << " --> New image estimate" 00896 << " | mean: " << mp_imageStatMean[0] 00897 << " | stdv: " << sqrt(mp_imageStatVariance[0]) 00898 << " | min/max: " << mp_imageStatMin[0] << "/" << mp_imageStatMax[0] << endl); 00899 } 00900 } 00901 } 00902 } 00903 00904 // End 00905 return 0; 00906 } 00907 00908 // ===================================================================== 00909 // --------------------------------------------------------------------- 00910 // --------------------------------------------------------------------- 00911 // ===================================================================== 00912 00913 FLTNB vOptimizer::ForwardProject(oProjectionLine* ap_Line, FLTNB* ap_image) 00914 { 00915 // Particular case for SPECT with attenuation correction 00916 if (m_dataType==TYPE_SPECT && m2p_attenuationImage[ap_Line->GetThreadNumber()]!=NULL) 00917 return ap_Line->ForwardProjectWithSPECTAttenuation( m2p_attenuationImage[ap_Line->GetThreadNumber()], ap_image ); 00918 // Otherwise call the function without attenuation 00919 else 00920 return ap_Line->ForwardProject( ap_image ); 00921 } 00922 00923 // ===================================================================== 00924 // --------------------------------------------------------------------- 00925 // --------------------------------------------------------------------- 00926 // ===================================================================== 00927 00928 void vOptimizer::BackwardProject(oProjectionLine* ap_Line, FLTNB* ap_image, FLTNB a_value) 00929 { 00930 // Particular case for SPECT with attenuation correction 00931 if (m_dataType==TYPE_SPECT && m2p_attenuationImage[ap_Line->GetThreadNumber()]!=NULL) 00932 ap_Line->BackwardProjectWithSPECTAttenuation( m2p_attenuationImage[ap_Line->GetThreadNumber()], ap_image, a_value ); 00933 // Otherwise call the function without attenuation 00934 else 00935 ap_Line->BackwardProject( ap_image, a_value ); 00936 } 00937 00938 // ===================================================================== 00939 // --------------------------------------------------------------------- 00940 // --------------------------------------------------------------------- 00941 // ===================================================================== 00942 00943 FLTNB vOptimizer::ComputeSensitivity( FLTNB**** a4p_sensitivityMatrix, 00944 int a_timeBasisFunction, int a_respBasisFunction, int a_cardBasisFunction, 00945 int a_voxel, int a_nbSubsets ) 00946 { 00947 // The sensitivity to be computed 00948 FLTNB sensitivity = 0.; 00949 // Loop over time frames 00950 for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++) 00951 { 00952 // Get frame/basis coefficient and continue if null 00953 FLTNB time_basis_coef = mp_ImageDimensionsAndQuantification->GetTimeBasisCoefficient(a_timeBasisFunction,fr); 00954 if (time_basis_coef==0.) continue; 00955 // Loop over respiratory gates 00956 for (int rg=0; rg<mp_ImageDimensionsAndQuantification->GetNbRespGates(); rg++) 00957 { 00958 // Get resp_gate/basis coefficient and continue if null 00959 FLTNB resp_basis_coef = mp_ImageDimensionsAndQuantification->GetRespBasisCoefficient(a_respBasisFunction,rg); 00960 if (resp_basis_coef==0.) continue; 00961 // Loop over cardiac gates 00962 for (int cg=0; cg<mp_ImageDimensionsAndQuantification->GetNbCardGates(); cg++) 00963 { 00964 // Get card_gate_basis coefficient and continue if null 00965 FLTNB card_basis_coef = mp_ImageDimensionsAndQuantification->GetCardBasisCoefficient(a_cardBasisFunction,cg); 00966 if (card_basis_coef==0.) continue; 00967 // Add actual weight contribution 00968 sensitivity += a4p_sensitivityMatrix[fr][rg][cg][a_voxel] * 00969 time_basis_coef * resp_basis_coef * card_basis_coef; 00970 } 00971 } 00972 } 00973 // If listmode, then divide by the number of subsets 00974 if (m_dataMode==MODE_LIST) sensitivity /= ((FLTNB)a_nbSubsets); 00975 // Return sensitivity 00976 return sensitivity; 00977 } 00978 00979 // ===================================================================== 00980 // --------------------------------------------------------------------- 00981 // --------------------------------------------------------------------- 00982 // =====================================================================
1.7.6.1
