128 for (
int bed=0 ; bed<
m_nbBeds ; bed++)
131 #pragma omp parallel for private(index) schedule(static,1) 137 th = omp_get_thread_num();
160 if (
m_verbose>=2)
Cout(
"iRCPGSAlgorithm::InitSpecificOptions() "<<a_specificOptions<<endl);
168 size_t colon = a_specificOptions.find_first_of(
":");
169 size_t comma = a_specificOptions.find_first_of(
",");
172 string list_options =
"";
173 string file_options =
"";
176 if (colon!=string::npos)
179 name = a_specificOptions.substr(0,colon);
181 file_options = a_specificOptions.substr(colon+1);
186 else if (comma!=string::npos)
189 name = a_specificOptions.substr(0,comma);
191 list_options = a_specificOptions.substr(comma+1);
199 name = a_specificOptions;
204 file_options = p_output_manager->
GetPathToConfigDir() +
"/algorithm/" + name +
".conf";
210 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> The given algorithm name is not RCPGS !" << endl);
215 if(!list_options.empty())
217 vector<string> option_name;
218 vector<string> option_value;
222 while ((pos_comma=list_options.find_first_of(
","))!=string::npos)
225 string sub_buf = list_options.substr(0,pos_comma);
227 size_t pos_equal = sub_buf.find_first_of(
"=");
228 if (pos_equal==string::npos || pos_equal==0)
230 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Syntax problem in algorithm parameters !" << endl);
233 option_name.push_back(sub_buf.substr(0,pos_equal));
234 option_value.push_back(sub_buf.substr(pos_equal+1));
235 list_options = list_options.substr(pos_comma+1);
239 if (list_options.length()>0)
241 size_t pos_equal = list_options.find_first_of(
"=");
242 if (pos_equal==string::npos || pos_equal==0)
244 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Syntax problem in algorithm parameters !" << endl);
247 option_name.push_back(list_options.substr(0,pos_equal));
248 option_value.push_back(list_options.substr(pos_equal+1));
252 vector<HPFLTNB> temp_multiModalNoiseSigma;
255 for(
size_t v=0;v<option_name.size();v++)
257 if (option_name.at(v)==
"ddCRP")
263 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Failed to read the list of options correctly!" << endl);
268 else if (option_name.at(v)==
"alpha")
274 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Failed to read the list of options correctly!" << endl);
279 else if (option_name.at(v)==
"gammaShape")
285 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Failed to read the list of options correctly!" << endl);
290 else if (option_name.at(v)==
"gammaRate")
296 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Failed to read the list of options correctly!" << endl);
301 else if (option_name.at(v)==
"backprojection")
307 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Failed to read the list of options correctly!" << endl);
312 else if (option_name.at(v)==
"multiModalNoiseSigma")
316 if (
ReadStringOption(option_value.at(v), option, 1,
",",
"multiModalNoiseSigma"))
318 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Failed to read the list of options correctly!" << endl);
321 temp_multiModalNoiseSigma.push_back(option[0]);
323 else if (option_name.at(v)==
"multiModalLag")
329 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Failed to read the list of options correctly!" << endl);
334 else if (option_name.at(v)==
"meanClusterVolumeMin")
338 if (
ReadStringOption(option_value.at(v), option, 1,
",",
"meanClusterVolumeMin"))
340 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Failed to read the list of options correctly!" << endl);
345 else if (option_name.at(v)==
"meanClusterVolumeMax")
349 if (
ReadStringOption(option_value.at(v), option, 1,
",",
"meanClusterVolumeMax"))
351 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Failed to read the list of options correctly!" << endl);
356 else if (option_name.at(v)==
"alphaIncrement")
362 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Failed to read the list of options correctly!" << endl);
370 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> The number of provided noise standard deviations does not equal the number of provided multimodal images! " << endl);
380 Cout(
"iRCPGSAlgorithm::InitSpecificOptions() -> Algorithm options read from command line "<<endl);
384 else if(!file_options.empty())
389 Cout(
"iRCPGSAlgorithm::InitSpecificOptions() -> Algorithm options read from the configuration file "<<endl);
396 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> The ddCRP alpha parameter must be non negative when ddCRP is used! " << endl);
399 if (m_ddCRP<0 || m_ddCRP>2)
401 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> The ddCRP type is not properly specified! " << endl);
404 if (m_backprojection<0 || m_backprojection>2)
406 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> The backprojection type is not properly specified! " << endl);
411 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Gamma prior parameters are not positive! " << endl);
418 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Wrong noise standard deviations for multimodal images! " << endl);
427 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Min/Max criteria for average cluster volume inconsistent! " << endl);
432 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Multiplicative update for ddCRP alpha must be >1. if min/max criteria for average cluster volume specified ! " << endl);
437 Cerr(
"***** iRCPGSAlgorithm::InitSpecificOptions() -> Options for adaptative alpha not compatible with alpha=0. ! " << endl);
449 Cout(
" as alpha=0, entering special mode when sampling the posterior clustering"<<endl);
456 else Cout(
"iRCPGSAlgorithm::InitSpecificOptions() -> ddCRP not used ! "<<endl);
457 Cout(
"iRCPGSAlgorithm::InitSpecificOptions() -> Backprojection type = "<<
m_backprojection<<endl);
465 string key_word =
"";
467 key_word =
"ddCRP type";
470 Cerr(
"***** iRCPGSAlgorithm::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word <<
"' keyword !" << endl);
474 key_word =
"ddCRP alpha";
477 Cerr(
"***** iRCPGSAlgorithm::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word <<
"' keyword !" << endl);
481 key_word =
"gamma prior shape";
484 Cerr(
"***** iRCPGSAlgorithm::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word <<
"' keyword !" << endl);
488 key_word =
"gamma prior rate";
491 Cerr(
"***** iRCPGSAlgorithm::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word <<
"' keyword !" << endl);
495 key_word =
"backprojection type";
498 Cerr(
"***** iRCPGSAlgorithm::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word <<
"' keyword !" << endl);
502 key_word =
"multimodal noise sigma";
505 Cerr(
"***** iRCPGSAlgorithm::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word <<
"' keyword !" << endl);
509 key_word =
"multimodal lag";
512 Cerr(
"***** iRCPGSAlgorithm::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word <<
"' keyword !" << endl);
516 key_word =
"mean cluster volume min";
519 Cerr(
"***** iRCPGSAlgorithm::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word <<
"' keyword !" << endl);
523 key_word =
"mean cluster volume max";
526 Cerr(
"***** iRCPGSAlgorithm::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word <<
"' keyword !" << endl);
530 key_word =
"ddCRP alpha increment";
533 Cerr(
"***** iRCPGSAlgorithm::ReadAndCheckConfigurationFile() -> Failed to get the '" << key_word <<
"' keyword !" << endl);
545 if (
m_verbose>=2)
Cout(
"iRCPGSAlgorithm::InitializeHelperVar() "<<endl);
596 for (
int bed=0 ; bed<
m_nbBeds ; bed++)
602 int64_t index_start = 0;
603 int64_t index_stop = 0;
604 m2p_DataFile[bed]->GetEventIndexStartAndStop(&index_start, &index_stop, 0, 1);
608 #pragma omp parallel for private(index) schedule(static, 1) 609 for (index=0 ; index<m2p_DataFile[bed]->GetSize() ; index++)
614 th = omp_get_thread_num();
617 vEvent*
event = m2p_DataFile[bed]->GetEvent(index, th);
626 for (
int e=0; e<(
INTNB)round(event->GetEventValue(b)); e++)
637 #pragma omp parallel for private(v) schedule(guided) 638 for (v=0; v<nbVoxels; v++)
690 if (
m_verbose>=2)
Cout(
"iRCPGSAlgorithm::ProcessMultiModalInfo() "<<endl);
719 Cerr(
"***** iRCPGSAlgorithm::StepBeforeIterationLoop() -> A problem occurred while calling StepBeforeIterationLoop() function !" << endl);
726 Cerr(
"***** vAlgorithm::StepBeforeIterationLoop() -> An error occurred while reading the initialization image !" << endl);
743 if (
m_verbose>=2)
Cout(
"iRCPGSAlgorithm::StepBeforeIterationLoop() --> Dividing the provided sensitivity image by the number of subsets "<<endl);
762 if (
m_verbose>=3)
Cout(
"iRCPGSAlgorithm::StepPreProcessInsideSubsetLoop ... " << endl);
794 if (
m_verbose>=3)
Cout(
"iRCPGSAlgorithm::StepPostProcessInsideSubsetLoop ... " << endl);
809 #pragma omp parallel for private(v) schedule(guided) 834 #pragma omp parallel for private(s) schedule(guided) 858 bool alpha_changed =
false;
863 alpha_changed =
true;
869 alpha_changed =
true;
875 Cout(
"iRCPGSAlgorithm::StepPostProcessInsideSubsetLoop() -> Updated ddCRP alpha to "<<
m_ddcrpAlpha<<endl);
886 stringstream temp_it; temp_it << a_iteration + 1;
887 stringstream temp_ss; temp_ss << a_subset + 1;
888 temp_sens.append(
"_it").append(temp_it.str()).append(
"_ss").append(temp_ss.str()).append(
"_sensitivity");
897 if (
m_verbose>=1)
Cout(
"iRCPGSAlgorithm::StepPostProcessInsideSubsetLoop() -> Save image at iteration " << a_iteration+1 <<
" for subset " << a_subset+1 << endl);
901 Cerr(
"***** iRCPGSAlgorithm::StepPostProcessInsideSubsetLoop() -> A problem occurred while saving images at iteration " << a_iteration+1 <<
" !" << endl);
918 if (
m_nbBeds>1)
Cout(
"iRCPGSAlgorithm::SampleConditionalCompleteData() -> Start loop over events for bed " << a_bed+1 << endl << flush);
919 else Cout(
"iRCPGSAlgorithm::SampleConditionalCompleteData() -> Start loop over events" << endl << flush);
931 cout <<
"0 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%" << endl;
932 cout <<
"|----|----|----|----|----|----|----|----|----|----|" << endl;
933 cout <<
"|" << flush;
936 int progression_percentage_old = 0;
937 int progression_nb_bars = 0;
938 uint64_t progression_printing_index = 0;
941 int64_t index_start = 0;
942 int64_t index_stop = 0;
948 MPI_Barrier(MPI_COMM_WORLD);
957 bool problem =
false;
971 #pragma omp parallel for private(index_temp) schedule(static, 1) 972 for ( index_temp = index_start ; index_temp < index_stop ; index_temp +=
mp_nbSubsets[a_iteration])
977 th = omp_get_thread_num();
981 INTNB index = index_temp;
987 if (progression_printing_index%1000==0)
989 int progression_percentage_new = ((int)( (((
float)(index-index_start+1))/((float)(index_stop-index_start)) ) * 100.));
990 if (progression_percentage_new>=progression_percentage_old+2)
992 int nb_steps = (progression_percentage_new-progression_percentage_old)/2;
993 for (
int i=0; i<nb_steps; i++)
995 cout <<
"-" << flush;
996 progression_nb_bars++;
998 progression_percentage_old += nb_steps*2;
1001 progression_printing_index++;
1011 Cerr(
"***** iRCPGSAlgorithm::SampleConditionalCompleteData() -> An error occured while getting the event from index " 1012 << index <<
" (thread " << th <<
") !" << endl);
1025 Cout(
"iRCPGSAlgorithm::SampleConditionalCompleteData() -> Step2: Check for Dynamic event (frame/gate switch, image-based deformation " << endl);
1028 int dynamic_switch_value =
mp_ID->
DynamicSwitch(index, event->GetTimeInMs(), a_bed, th);
1039 bool emptyEvent =
true;
1040 for (
int b=0; b<
event->GetNbValueBins(); b++)
1041 if (event->GetEventValue(b)>0.)
1047 if (emptyEvent)
continue;
1052 if (
m_verbose>=4)
Cout(
"iRCPGSAlgorithm::SampleConditionalCompleteData() -> Step3: Compute the projection line " << endl);
1057 Cerr(
"***** iRCPGSAlgorithm::SampleConditionalCompleteData() -> A problem occured while computing the projection line !" << endl);
1072 assert(timeFrame==0);assert(respGate==0);assert(cardGate==0);
1079 if (time_basis_coef==0.)
continue;
1085 if (resp_basis_coef==0.)
continue;
1091 if (card_basis_coef==0.)
continue;
1093 FLTNB global_basis_coef = time_basis_coef * resp_basis_coef * card_basis_coef;
1094 assert(global_basis_coef==1.);
1102 FLTNB temp_eventValue =
event->GetEventValue(b);
1103 if (temp_eventValue>0.)
1108 INTNB numCounts = (
INTNB)round(temp_eventValue);
1111 assert(round(temp_eventValue)==temp_eventValue);
1114 for (
int e=0; e<numCounts; e++)
1135 mp_clusterN[th][previous_voxel_current_cluster] -= 1.;
1136 assert(
mp_clusterN[th][previous_voxel_current_cluster]>-0.001);
1145 vector<HPFLTNB> voxelProbabilities(current_nb_voxels+1,0.);
1146 HPFLTNB probabilities_sum = 0.;
1148 INTNB temp_voxelIndex = -1;
1153 for (
int vl=0; vl<current_nb_voxels; vl++)
1158 probabilities_sum += temp_value;
1159 voxelProbabilities.at(vl+1) = temp_value;
1165 for (
int vl=0; vl<current_nb_voxels; vl++)
1180 probabilities_sum += temp_value;
1181 voxelProbabilities.at(vl) = temp_value;
1188 voxelProbabilities[current_nb_voxels] = temp_value;
1189 probabilities_sum += temp_value;
1192 assert(!(probabilities_sum<0.));
1194 if (probabilities_sum>0.)
1201 INTNB sampledIndex = 0;
1202 for (sampledIndex=current_nb_voxels; sampledIndex>=0; sampledIndex--)
1204 cumulative_sum += voxelProbabilities[sampledIndex];
1205 if (cumulative_sum>current_random)
break;
1209 if (sampledIndex>=0 && sampledIndex<current_nb_voxels)
1244 vector<HPFLTNB> voxelProbabilities(current_nb_voxels+1,0.);
1245 HPFLTNB probabilities_sum = 0.;
1248 for (
int vl=0; vl<current_nb_voxels; vl++)
1255 probabilities_sum += temp_value;
1256 voxelProbabilities.at(vl) = temp_value;
1261 voxelProbabilities[current_nb_voxels] = temp_value;
1262 probabilities_sum += temp_value;
1264 assert(!(probabilities_sum<0.));
1268 if (probabilities_sum>0.)
1271 INTNB numCounts = (
INTNB)round(temp_eventValue);
1274 if (fabs(round(temp_eventValue)-temp_eventValue)>1.e-7)
1278 numCounts = (
INTNB)ceil(temp_eventValue);
1282 numCounts = (
INTNB)floor(temp_eventValue);
1288 for (
int e=0; e<numCounts; e++)
1295 INTNB sampledIndex = 0;
1297 for (sampledIndex=current_nb_voxels; sampledIndex>=0; sampledIndex--)
1299 cumulative_sum += voxelProbabilities[sampledIndex];
1300 if (cumulative_sum>current_random)
break;
1304 if (sampledIndex>=0 && sampledIndex<current_nb_voxels)
1308 else if (sampledIndex<0)
1311 Cout(
"Warning: Nothing sampled, sampledIndex "<<sampledIndex<<
", current_nb_voxels "<<current_nb_voxels<<endl);
1327 MPI_Barrier(MPI_COMM_WORLD);
1338 int progression_total_bars = 49;
1339 for (
int i=0; i<progression_total_bars-progression_nb_bars; i++) cout <<
"-";
1340 cout <<
"|" << endl;
1346 Cerr(
"***** iRCPGSAlgorithm::SampleConditionalCompleteData() -> A problem occured inside the parallel loop over events !" << endl);
1358 if (
m_verbose>=2)
Cout(
"iRCPGSAlgorithm::ComputeSumsPerClusters()... " << endl);
1372 #pragma omp parallel for private(vp) schedule(guided) 1378 if (multimodal_info)
1404 if (multimodal_info)
1423 if (
m_verbose>=3)
Cout(
"iRCPGSAlgorithm::UpdateVisitedVoxels() -> Tick visited voxels based on sensitivity" << endl);
1443 if (
m_verbose>=2)
Cout(
"iRCPGSAlgorithm::SampleConditionalClustering()... " << endl);
1451 vector<INTNB> fellow_voxels;
1456 vector<HPFLTNB> probabilities;
1461 vector<INTNB> voxels_after_cut;
1471 for(
int v=0; v<dim_3D; v++)
1483 assert( temp_new_cluster<=dim_3D );
1484 assert(temp_new_cluster!=previous_cluster);
1492 voxels_after_cut.resize(0);
1495 function<void (int)> explorePreviousLinks= [
this,&explorePreviousLinks,current_voxel,&voxels_after_cut,&loop] (
int index)
1497 voxels_after_cut.push_back(index);
1500 if (ip!=current_voxel) explorePreviousLinks(ip);
1504 explorePreviousLinks(current_voxel);
1505 assert( voxels_after_cut.size()>0 );
1508 size_t previous_next_link =
mp_nextLink[current_voxel];
1522 INTNB temp_count = 0;
1525 for(
int nv : voxels_after_cut)
1531 if (multimodal_info)
1542 assert (temp_count_N>=0. && temp_count_sens>=0.);
1543 assert (temp_count>=0);
1561 if (multimodal_info)
1577 probabilities.resize(0);
1582 for(
size_t fv=1; fv<fellow_voxels.size(); fv++)
1585 if(fellow_cluster!=temp_new_cluster)
1593 assert (!(n0<0. || n1<0. || s0<0. || s1<0.));
1605 assert (c0>0. && c1>0.);
1606 HPFLTNB help1 = c0*c1/(c0+c1);
1612 assert (a0>=0. && a1>=0.);
1625 assert (!std::isnan(temp));
1626 probabilities.push_back(temp);
1637 probabilities.push_back (0.);
1645 for(
size_t p=0; p<probabilities.size(); p++)
1647 probabilities[p] = exp(probabilities[p]);
1648 assert (!std::isnan(probabilities[p]));
1649 assert (!std::isinf(probabilities[p]));
1650 probabilities_sum += probabilities[p];
1651 assert (!std::isinf(probabilities_sum));
1658 if (probabilities_sum>std::numeric_limits<HPFLTNB>::min())
1662 for(p=0; p<(
INTNB)probabilities.size(); p++)
1664 cumulative_sum += probabilities[p];
1665 if(cumulative_sum>current_random)
break;
1677 assert(p<(
INTNB)probabilities.size());
1684 int new_next_link = fellow_voxels.at(p);
1695 if (new_next_link_cluster!=temp_new_cluster)
1701 assert(
mp_clusterN[0][new_next_link_cluster]>=0.);
1704 if (multimodal_info)
1718 for(
size_t fv=0; fv<voxels_after_cut.size(); fv++)
1740 if (
m_verbose>=2)
Cout(
"iRCPGSAlgorithm::SampleConditionalClusterIntensity()... " << endl);
1745 #pragma omp parallel for private(vp) schedule(guided) 1755 nb_relevant_voxels ++;
1759 if (
m_ddCRP==0) assert(index==v);
1780 Cout(
" --> Mean number of voxels per cluster = "<<((
FLTNB)nb_relevant_voxels/(
FLTNB)nb_clusters)<<endl);
1800 a_fellow_voxels.resize(0);
1802 a_fellow_voxels.push_back(a_current_voxel);
1807 const int ix = a_current_voxel%dim_X;
1813 if (y_low) a_fellow_voxels.push_back(a_current_voxel-dim_X);
1814 if (y_high) a_fellow_voxels.push_back(a_current_voxel+dim_X);
1815 if (x_low) a_fellow_voxels.push_back(a_current_voxel-1);
1816 if (x_high) a_fellow_voxels.push_back(a_current_voxel+1);
1818 for(
size_t i=0;i<a_fellow_voxels.size();i++) assert(a_fellow_voxels.at(i)>=0 && a_fellow_voxels.at(i)<dim_3D );
1819 assert(a_fellow_voxels.size()>0);
1824 const int current_voxel_xy = a_current_voxel%(dim_XY);
1826 const int ix = current_voxel_xy%dim_X;
1835 if (y_low) a_fellow_voxels.push_back(a_current_voxel-dim_X);
1836 if (y_high) a_fellow_voxels.push_back(a_current_voxel+dim_X);
1837 if (x_low) a_fellow_voxels.push_back(a_current_voxel-1);
1838 if (x_high) a_fellow_voxels.push_back(a_current_voxel+1);
1839 if (z_low) a_fellow_voxels.push_back(a_current_voxel-dim_XY);
1840 if (z_high) a_fellow_voxels.push_back(a_current_voxel+dim_XY);
1842 for(
size_t i=0;i<a_fellow_voxels.size();i++) assert(a_fellow_voxels.at(i)>=0 && a_fellow_voxels.at(i)<dim_3D );
1843 assert(a_fellow_voxels.size()>0);
1856 #pragma omp parallel for private(vp) schedule(guided) 1883 Cerr(
"***** iRCPGSAlgorithm::StepAfterIterationLoop() -> A problem occurred while calling StepAfterIterationLoop() function !" << endl);
1886 if (
m_verbose>=2)
Cout(
"iRCPGSAlgorithm::StepAfterIterationLoop ... " << endl);
1902 if (
m_verbose>=3)
Cout(
"iRCPGSAlgorithm::StepBeforeSubsetLoop ... " << endl);
1912 if (
m_verbose>=3)
Cout(
"iRCPGSAlgorithm::StepAfterSubsetLoop() -> Clean never visited voxels and save images if needed" << endl);
1919 if (
m_verbose>=1)
Cout(
"iRCPGSAlgorithm::StepAfterSubsetLoop() -> Save image at iteration " << a_iteration+1 << endl);
1923 Cerr(
"***** iRCPGSAlgorithm::StepAfterSubsetLoop() -> A problem occurred while saving images at iteration " << a_iteration+1 <<
" !" << endl);
1949 Cerr(
"***** iRCPGSAlgorithm::StepImageOutput() -> A problem occurred while applying output FOV masking !" << endl);
1955 Cerr(
"***** iRCPGSAlgorithm::StepImageOutput() -> A problem occurred while applying output flip !" << endl);
1966 Cerr(
"***** iRCPGSAlgorithm::StepImageOutput() -> A problem occurred while saving output image !" << endl);
1970 if (((a_iteration+1)%50)==0)
1981 if (a_iteration >= 0)
1983 stringstream ss; ss << a_iteration + 1;
1984 data_file.append(
"_it").append(ss.str());
1987 data_file.append(
".img");
1994 FILE* fout = fopen(data_file.c_str(),
"wb");
1997 Cerr(
"***** iRCPGSAlgorithm::StepImageOutput() -> Failed to create output file '" << data_file <<
"' !" << endl);
2008 nb_data += fwrite(&buffer,
sizeof(
INTNB),1,fout);
2017 Cerr(
"***** iRCPGSAlgorithm::StepImageOutput() -> Failed to write all data into the output file '" << data_file <<
"' !" << endl);
2034 cout <<
"Random Clustering Prior - Gibbs Sampler : probabilistic Bayesian inference algorithm (no optimization)" << endl;
2035 cout <<
"M. Filipovic et al : PET reconstruction of the posterior image probability, including multimodal images." << endl;
2037 cout <<
"The algorithm is launched using -prob, either with command line options (-prob RCPGS:option1name=option1value,option2name=option2value,...)" << endl;
2038 cout <<
"or with a config file (-prob RCPGS:config_file.conf), see also the default config file for examples of option values." << endl;
2040 cout <<
"Required options :" << endl;
2041 cout <<
" ddCRP = type of ddCRP prior (0 = no ddCRP, 1 = original ddCRP (recommended), 2 = modified ddCRP)" << endl;
2042 cout <<
" alpha = ddCRP hyperparameter, the unnormalized probability of drawing a self link" << endl;
2043 cout <<
" gammaShape = Gamma prior distribution shape parameter (0.5 recommended)" << endl;
2044 cout <<
" gammaRate = Gamma prior distribution rate parameter (1.e-18 recommended)" << endl;
2045 cout <<
" backprojection = type of multinomial backprojection of the current iteration/subset data (0 = the previous backprojection state is cleared, as for ML-EM (recommended), 1 = update of the previous backprojection state, 2 = backprojection marginalized over cluster intensity, implies update of the previous backprojection state )" << endl;
2046 cout <<
" multiModalNoiseSigma = standard deviation of Gaussian noise in the provided additional image from another modality, must be repeated as many times as the -multimodal option" << endl;
2047 cout <<
" multiModalLag = number of iterations after which the multimodal images start affecting voxels clustering" << endl;
2049 cout <<
"Optional options (if one of these is specified, all the others must be specified as well):" << endl;
2050 cout <<
" meanClusterVolumeMin = minimum threshold for average cluster volume (in mm3), used for tuning ddCRP alpha automatically through iterations" << endl;
2051 cout <<
" meanClusterVolumeMax = maximum threshold for average cluster volume (in mm3), used for tuning ddCRP alpha automatically through iterations" << endl;
2052 cout <<
" alphaIncrement = multiplicative increment for tuning ddCRP alpha through iterations" << endl;
void StopIterativeDataUpdateStep1(int a_thread)
Stop the timer for duration of iterative data update step 1.
FLTNB **** m4p_forwardImage
int UpdateVisitedVoxels()
Check for voxels which do not contribute to recorded observed data.
int GetCurrentCardImage(int a_th)
call the eponym function from the oDynamicDataManager object
int SampleConditionalClustering(int a_iteration)
Gibbs sampler : second conditional probability (ddCRP links)
HPFLTNB m_ddcrpAlphaIncrement
void GetEventIndexStartAndStop(int64_t *ap_indexStart, int64_t *ap_indexStop, int a_subsetNum=0, int a_NbSubsets=1)
Compute the index start and stop of the events loop with respect to the current subset and MPI size a...
static sRandomNumberGenerator * GetInstance()
Instanciate the singleton object and Initialize member variables if not already done, return a pointer to this object otherwise.
void DeallocateBackwardImageFromDynamicBasis()
Free memory for the backward image matrices.
int InitSpecificOptions(string a_specificOptions)
int GetNbCardBasisFunctions()
Get the number of cardiac basis functions.
int ComputeFellowVoxelsList(vector< INTNB > &a_fellow_voxels, int a_current_voxel)
Fill the input list of fellow voxels for the current type of neighbourhood and the input voxel...
FLTNB GetVoxSizeX()
Get the voxel's size along the X axis, in mm.
#define DYNAMIC_SWITCH_CONTINUE
HPFLTNB ** mpp_clusterMultiModal
FLTNB * mp_permanentBackwardImage
int StepPreProcessInsideSubsetLoop(int a_iteration, int a_subset)
This function is called right after starting the data subset loop.
HPFLTNB m_currentMeanClusterVolume
void Reduce()
Merge parallel results into the matrix of the backward image matrix of the first thread. Also for MPI.
FLTNB GetVoxSizeZ()
Get the voxel's size along the Z axis, in mm.
void StartCustomStep(int a_thread, int a_step)
Start the timer for duration of custom step of the given index for the given thread.
FLTNB GetRespBasisCoefficient(int a_respBasisFunction, int a_respGate)
Get the respiratory basis coefficients for the provided respiratory gate and basis function...
bool IsLoadedSensitivity()
int DynamicSwitch(int64_t a_currentEventIndex, uint32_t a_currentTime, int a_bed, int a_th)
Call the eponym function from the oDynamicDataManager class using the member object.
HPFLTNB * mp_clusterValues
int ApplyOutputFOVMasking()
Mask the outside of the transaxial FOV based on the m_fovOutPercent.
int ReadConfigurationFile(const string &a_configurationFile)
oProjectorManager * mp_ProjectorManager
int StepInnerLoopInsideSubsetLoop(int a_iteration, int a_subset, int a_bed)
This function is called inside the subset loop. It contains the core operations of the algorithm an...
bool m_saveSensitivityHistoFlag
void StopIterativeDataUpdateStep3(int a_thread)
Stop the timer for duration of iterative data update step 3.
void InstantiateBackwardImageFromDynamicBasis(int a_nbBackwardImages)
Allocate memory for the backward image matrices and set the number of backward images for the whole c...
int GetNbMultiModalImages()
Get the number of additional multimodal images.
HPFLTNB GenerateExtraRdmNber(int a_nb=0)
Generate a random number using the specified additional not thread safe random generator, for use in sequential parts of an otherwise multithreaded code.
INTNB ** mp_listEventsIndices
int ApplyMaskToSensitivity()
Apply the mask to the sensitivity image (only for the first thread, the image must be reduced beforeh...
int GetNbTimeBasisFunctions()
Get the number of time basis functions.
HPFLTNB * mp_multiModalParam
static sOutputManager * GetInstance()
Instanciate the singleton object and Initialize member variables if not already done, return a pointer to this object otherwise.
int InitImage(const string &a_pathToInitialImage, FLTNB a_value)
Initialize the main image, either using:
int ComputeSumsPerClusters(int a_iteration)
Compute sums of voxel features for each cluster.
void StartIterativeDataUpdateStep1(int a_thread)
Start the timer for duration of iterative data update step 1.
oImageSpace * mp_ImageSpace
void StartIterativeDataUpdateStep2(int a_thread)
Start the timer for duration of iterative data update step 2.
FLTNB * mp_visitedVoxelsImage
vEvent * GetEvent(int64_t a_eventIndex, int a_th=0)
vDataFile ** m2p_DataFile
INTNB GetVoxelIndex(int a_direction, int a_TOFBin, INTNB a_voxelInLine)
This function is used to get the contributing voxel index of the provided direction, TOF bin and voxel rank.
bool * mp_listRelevantVoxelIndices
void PrepareForwardImage()
Copy current image matrix in the forward-image buffer matrix.
FLTNB GetVoxSizeY()
Get the voxel's size along the Y axis, in mm.
HPFLTNB m_meanClusterVolumeMax
const string & GetPathName()
FLTNB ****** m6p_backwardImage
HPFLTNB m_meanClusterVolumeMin
void CleanNeverVisitedVoxels()
Based on the visitedVoxelsImage, clean the never visited voxels in the image. This function must be c...
Singleton class that manages output writing on disk (images, sinograms, etc). It also manages loggi...
const string & GetPathToConfigDir()
Return the path to the CASTOR config directory.
string m_pathToInitialImg
FLTNB ** m2p_multiModalImage
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...
int StepPostProcessInsideSubsetLoop(int a_iteration, int a_subset)
int SaveOutputImage(int a_iteration, int a_subset=-1)
Call the interfile function to write output image on disk.
int StepBeforeIterationLoop()
This function is called at the beginning of the RunCPU function.
void ResetCurrentDynamicIndices()
Call the eponym function from the oDynamicDataManager class using the member object.
virtual int StepBeforeIterationLoop()
This function is called at the beginning of the RunCPU function.
bool m_saveImageAfterSubsets
INTNB GetCurrentNbVoxels(int a_direction, int a_TOFBin)
This function is used to get the current number of contributing voxels to the line.
void StartIterativeDataUpdateStep3(int a_thread)
Start the timer for duration of iterative data update step 3.
oProjectionLine * ComputeProjectionLine(vEvent *ap_Event, int a_th)
This function is used to compute system matrix elements from the associated projector or pre-computed...
void StopIterativeDataUpdateStep2(int a_thread)
Stop the timer for duration of iterative data update step 2.
FLTNB GetQuantificationFactor(int a_bed, int a_frame, int a_respGate, int a_cardGate)
Get the quantification factor corresponding to the provided bed, frame, respiratory and cardiac gates...
INTNB GetNbVoxXY()
Get the number of voxels in a slice.
bool IsLoadedMultiModal()
int SampleConditionalCompleteData(int a_iteration, int a_subset, int a_bed)
Gibbs sampler : first conditional probability (backprojection)
void ApplyBedOffset(int a_bed)
Compute the bed offset from the provided bed index and apply it to all projection lines...
HPFLTNB * mp_clusterSensitivity
int StepAfterSubsetLoop(int a_iteration)
This function is called after finishing the data subset loop.
HPFLTNB GenerateRdmNber()
Generate a random number for the thread which index is recovered from the OpenMP function.
FLTNB GetTimeBasisCoefficient(int a_timeBasisFunction, int a_timeFrame)
Get the time basis coefficients for the provided frame and basis function.
bool NotEmptyLine()
This function is used to know if the line contains any voxel contribution.
#define KEYWORD_MANDATORY
int StepBeforeSubsetLoop(int a_iteration)
This function is called before starting the data subset loop.
int InitializeHelperVar()
Allocate and initialize temporary variables, after the main variables have been initialized Assumptio...
INTNB **** m4p_EventsBackprojectionTrace
This class is designed to manage and store system matrix elements associated to a vEvent...
const string & GetBaseName()
int ApplyOutputFlip()
Just flip the output image.
oImageDimensionsAndQuantification * mp_ID
virtual int StepAfterIterationLoop()
This function is called at the end of the RunCPU function.
HPFLTNB * temp_count_multimodal
This is the base class for reconstructions, containing a framework with iteration and data subset loo...
HPFLTNB * mp_multiModalNoiseSigma
Mother class for the Event objects.
int GetNbThreadsForImageComputation()
Get the number of threads used for image operations.
INTNB GetNbVoxXYZ()
Get the total number of voxels.
vector< INTNB > * mpv_parentLinks
void ShowHelpSpecific()
Show help for the child algorithm.
int GetNbTOFBins()
This function is used to get the number of TOF bins in use.
FLTNB * AllocateMiscellaneousImage()
Allocate a new miscellaneous image on m2p_miscellaneousImages and return the pointer to this image...
int GetNbThreadsForProjection()
Get the number of threads used for projections.
FLTNB GetVoxelWeights(int a_direction, int a_TOFBin, INTNB a_voxelInLine)
This function is used to get the contributing voxel weight of the provided direction, TOF bin and voxel rank.
vector< INTNB > mv_newClusters
INTNB GetNbVoxX()
Get the number of voxels along the X axis.
bool * mp_outputIterations
INTNB * mp_voxelClusterMapping
int ProcessMultiModalInfo()
Check input multimodal images.
INTNB GetNbVoxZ()
Get the number of voxels along the Z axis.
int SampleConditionalClusterIntensity()
Gibbs sampler : third conditional probability (cluster intensity)
void StopCustomStep(int a_thread, int a_step)
Stop the timer for duration of custom step of the given index for the given thread.
This class is designed to manage some profiling of the code.
FLTNB GetFrameDurationInSec(int a_bed, int a_frame)
Get the frame duration for the given bed, in seconds as a FLTNB.
int StepAfterIterationLoop()
This function is called at the end of the RunCPU function.
int SaveSensitivityImage(const string &a_pathToSensitivityImage)
Call the interfile function to write the sensitivity image on disk.
void InitBackwardImage()
Initialize each voxel of the backward images to 0, also for sensitivity if not loaded (estimated on t...
INTNB * mp_listVoxelIndices
int GetCurrentRespImage(int a_th)
call the eponym function from the oDynamicDataManager object
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.
static sChronoManager * GetInstance()
Instantiate the singleton if not already done, then return the pointer to its instance.
FLTNB ***** m5p_sensitivity
int GetMPIRank()
Get the MPI instance number (the rank)
int GenerateCurrentImage()
Generate the current image estimation from the current sample of partition/clustering and cluster int...
int GetCurrentTimeFrame(int a_th)
call the eponym function from the oDynamicDataManager object
int StepImageOutput(int a_iteration, int a_subset=-1)
This function deals with everything about saving output images from the reconstruction.
void Display()
Display the results of the duration buffers.
FLTNB GetCardBasisCoefficient(int a_cardBasisFunction, int a_cardGate)
Get the cardiac basis coefficients for the provided cardiac gate and basis function.