From a32dac94273733f255bcd96f76a1f2bdc1545d77 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 15 May 2024 14:23:37 +0100 Subject: [PATCH] add TOF loop in computation of Hessian for projdata Fixes #1305 --- ...ihoodWithLinearModelForMeanAndProjData.cxx | 136 +++++++++++------- ...ihoodWithLinearModelForMeanAndProjData.cxx | 79 ++++++---- 2 files changed, 136 insertions(+), 79 deletions(-) diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx index 3cbd5832a..8c634e500 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx @@ -912,6 +912,30 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData::get_exam_info_up return exam_info_uptr; } +static std::vector +find_basic_viewgram_indices_in_subset(const ProjDataInfo& proj_data_info, + const DataSymmetriesForViewSegmentNumbers& symmetries, + const int min_segment_num, + const int max_segment_num, + const int subset_num, + const int num_subsets) +{ + const std::vector vs_nums_to_process = detail::find_basic_vs_nums_in_subset( + proj_data_info, symmetries, min_segment_num, max_segment_num, subset_num, num_subsets); + + std::vector vg_idx_to_process; + for (auto vs_num : vs_nums_to_process) + { + for (int k = proj_data_info.get_min_tof_pos_num(); k <= proj_data_info.get_max_tof_pos_num(); ++k) + { + ViewgramIndices viewgram_idx = vs_num; + viewgram_idx.timing_pos_num() = k; + vg_idx_to_process.push_back(viewgram_idx); + } + } + return vg_idx_to_process; +} + template Succeeded PoissonLogLikelihoodWithLinearModelForMeanAndProjData< @@ -938,35 +962,36 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData< this->get_projector_pair().get_forward_projector_sptr()->set_input(input); this->get_projector_pair().get_back_projector_sptr()->start_accumulating_in_new_target(); - const std::vector vs_nums_to_process - = detail::find_basic_vs_nums_in_subset(*this->get_proj_data().get_proj_data_info_sptr(), - *symmetries_sptr, - -this->get_max_segment_num_to_process(), - this->get_max_segment_num_to_process(), - subset_num, - this->get_num_subsets()); + const std::vector vg_idx_to_process + = find_basic_viewgram_indices_in_subset(*this->get_proj_data().get_proj_data_info_sptr(), + *symmetries_sptr, + -this->get_max_segment_num_to_process(), + this->get_max_segment_num_to_process(), + subset_num, + this->get_num_subsets()); info("Forward projecting input image.", 2); #ifdef STIR_OPENMP # pragma omp parallel for schedule(dynamic) #endif // note: older versions of openmp need an int as loop - for (int i = 0; i < static_cast(vs_nums_to_process.size()); ++i) + for (int i = 0; i < static_cast(vg_idx_to_process.size()); ++i) { + const auto viewgram_idx = vg_idx_to_process[i]; + { #ifdef STIR_OPENMP - const int thread_num = omp_get_thread_num(); - info(boost::format("Thread %d/%d calculating segment_num: %d, view_num: %d") % thread_num % omp_get_num_threads() - % vs_nums_to_process[i].segment_num() % vs_nums_to_process[i].view_num(), - 2); + const int thread_num = omp_get_thread_num(); + const int num_threads = omp_get_num_threads(); #else - info(boost::format("calculating segment_num: %d, view_num: %d") % vs_nums_to_process[i].segment_num() - % vs_nums_to_process[i].view_num(), - 2); + const int thread_num = 0; + const int num_threads = 1; #endif - const ViewSegmentNumbers view_segment_num = vs_nums_to_process[i]; - + info(boost::format("Thread %d/%d calculating segment_num: %d, view_num: %d, TOF: %d") % thread_num % num_threads + % viewgram_idx.segment_num() % viewgram_idx.view_num() % viewgram_idx.timing_pos_num(), + 2); + } // first compute data-term: y*norm^2 - RelatedViewgrams viewgrams = this->get_proj_data().get_related_viewgrams(view_segment_num, symmetries_sptr); + RelatedViewgrams viewgrams = this->get_proj_data().get_related_viewgrams(viewgram_idx, symmetries_sptr); // TODO add 1 for 1/(y+1) approximation this->get_normalisation().apply(viewgrams); @@ -978,7 +1003,7 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData< RelatedViewgrams tmp_viewgrams; // set tmp_viewgrams to geometric forward projection of input { - tmp_viewgrams = this->get_proj_data().get_empty_related_viewgrams(view_segment_num, symmetries_sptr); + tmp_viewgrams = this->get_proj_data().get_empty_related_viewgrams(viewgram_idx, symmetries_sptr); this->get_projector_pair().get_forward_projector_sptr()->forward_project(tmp_viewgrams); } @@ -1046,23 +1071,23 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData::actual_accumulat this->get_projector_pair().get_forward_projector_sptr()->set_input(input); this->get_projector_pair().get_back_projector_sptr()->start_accumulating_in_new_target(); - const std::vector vs_nums_to_process - = detail::find_basic_vs_nums_in_subset(*this->get_proj_data().get_proj_data_info_sptr(), - *symmetries_sptr, - -this->get_max_segment_num_to_process(), - this->get_max_segment_num_to_process(), - subset_num, - this->get_num_subsets()); + const std::vector vg_idx_to_process + = find_basic_viewgram_indices_in_subset(*this->get_proj_data().get_proj_data_info_sptr(), + *symmetries_sptr, + -this->get_max_segment_num_to_process(), + this->get_max_segment_num_to_process(), + subset_num, + this->get_num_subsets()); // Create and populate the input_viewgrams_vec with empty values. - // This is needed to make the order of the vector correct w.r.t vs_nums_to_process. + // This is needed to make the order of the vector correct w.r.t vg_idx_to_process. // OMP may mess this up - // Try: std::vector> input_viewgrams_vec(vs_nums_to_process.size()); + // Try: std::vector> input_viewgrams_vec(vg_idx_to_process.size()); std::vector> input_viewgrams_vec; - for (int i = 0; i < static_cast(vs_nums_to_process.size()); ++i) + for (int i = 0; i < static_cast(vg_idx_to_process.size()); ++i) { - const ViewSegmentNumbers view_segment_num = vs_nums_to_process[i]; - input_viewgrams_vec.push_back(this->get_proj_data().get_empty_related_viewgrams(view_segment_num, symmetries_sptr)); + const auto viewgram_idx = vg_idx_to_process[i]; + input_viewgrams_vec.push_back(this->get_proj_data().get_empty_related_viewgrams(viewgram_idx, symmetries_sptr)); } // Forward project input image @@ -1070,19 +1095,22 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData::actual_accumulat #ifdef STIR_OPENMP # pragma omp parallel for schedule(dynamic) #endif - for (int i = 0; i < static_cast(vs_nums_to_process.size()); ++i) - { // Loop over eah of the viewgrams in input_viewgrams_vec, forward projecting input into them + for (int i = 0; i < static_cast(vg_idx_to_process.size()); ++i) + { // Loop over each of the viewgrams in input_viewgrams_vec, forward projecting input into them + const auto viewgram_idx = vg_idx_to_process[i]; + { #ifdef STIR_OPENMP - const int thread_num = omp_get_thread_num(); - info(boost::format("Thread %d/%d calculating segment_num: %d, view_num: %d") % thread_num % omp_get_num_threads() - % vs_nums_to_process[i].segment_num() % vs_nums_to_process[i].view_num(), - 2); + const int thread_num = omp_get_thread_num(); + const int num_threads = omp_get_num_threads(); #else - info(boost::format("calculating segment_num: %d, view_num: %d") % vs_nums_to_process[i].segment_num() - % vs_nums_to_process[i].view_num(), - 2); + const int thread_num = 0; + const int num_threads = 1; #endif - input_viewgrams_vec[i] = this->get_proj_data().get_empty_related_viewgrams(vs_nums_to_process[i], symmetries_sptr); + info(boost::format("Thread %d/%d calculating segment_num: %d, view_num: %d, TOF: %d") % thread_num % num_threads + % viewgram_idx.segment_num() % viewgram_idx.view_num() % viewgram_idx.timing_pos_num(), + 2); + } + input_viewgrams_vec[i] = this->get_proj_data().get_empty_related_viewgrams(viewgram_idx, symmetries_sptr); this->get_projector_pair().get_forward_projector_sptr()->forward_project(input_viewgrams_vec[i]); } @@ -1091,35 +1119,37 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData::actual_accumulat #ifdef STIR_OPENMP # pragma omp parallel for schedule(dynamic) #endif - for (int i = 0; i < static_cast(vs_nums_to_process.size()); ++i) + for (int i = 0; i < static_cast(vg_idx_to_process.size()); ++i) { + const auto viewgram_idx = vg_idx_to_process[i]; + { #ifdef STIR_OPENMP - const int thread_num = omp_get_thread_num(); - info(boost::format("Thread %d/%d calculating segment_num: %d, view_num: %d") % thread_num % omp_get_num_threads() - % vs_nums_to_process[i].segment_num() % vs_nums_to_process[i].view_num(), - 2); + const int thread_num = omp_get_thread_num(); + const int num_threads = omp_get_num_threads(); #else - info(boost::format("calculating segment_num: %d, view_num: %d") % vs_nums_to_process[i].segment_num() - % vs_nums_to_process[i].view_num(), - 2); + const int thread_num = 0; + const int num_threads = 1; #endif + info(boost::format("Thread %d/%d calculating segment_num: %d, view_num: %d, TOF: %d") % thread_num % num_threads + % viewgram_idx.segment_num() % viewgram_idx.view_num() % viewgram_idx.timing_pos_num(), + 2); + } // Compute ybar_sq_viewgram = [ F(current_image_est) + additive ]^2 RelatedViewgrams ybar_sq_viewgram; { - ybar_sq_viewgram = this->get_proj_data().get_empty_related_viewgrams(vs_nums_to_process[i], symmetries_sptr); + ybar_sq_viewgram = this->get_proj_data().get_empty_related_viewgrams(vg_idx_to_process[i], symmetries_sptr); this->get_projector_pair().get_forward_projector_sptr()->forward_project(ybar_sq_viewgram); // add additive sinogram to forward projection if (!(is_null_ptr(this->get_additive_proj_data_sptr()))) - ybar_sq_viewgram += this->get_additive_proj_data().get_related_viewgrams(vs_nums_to_process[i], symmetries_sptr); + ybar_sq_viewgram += this->get_additive_proj_data().get_related_viewgrams(vg_idx_to_process[i], symmetries_sptr); // square ybar ybar_sq_viewgram *= ybar_sq_viewgram; } // Compute: final_viewgram * F(input) / ybar_sq_viewgram // final_viewgram starts as measured data - RelatedViewgrams final_viewgram - = this->get_proj_data().get_related_viewgrams(vs_nums_to_process[i], symmetries_sptr); + RelatedViewgrams final_viewgram = this->get_proj_data().get_related_viewgrams(vg_idx_to_process[i], symmetries_sptr); { // Mult input_viewgram final_viewgram *= input_viewgrams_vec[i]; diff --git a/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx b/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx index 73d47fece..45c111d0e 100644 --- a/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx +++ b/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx @@ -94,7 +94,7 @@ class PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests */ PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests(char const* const proj_data_filename = 0, char const* const density_filename = 0); - void construct_input_data(shared_ptr& density_sptr); + void construct_input_data(shared_ptr& density_sptr, const bool TOF_or_not); void run_tests() override; @@ -174,18 +174,37 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::test_approximate_Hes } void -PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::construct_input_data(shared_ptr& density_sptr) +PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::construct_input_data(shared_ptr& density_sptr, + const bool TOF_or_not) { if (this->proj_data_filename == 0) { + shared_ptr proj_data_info_sptr; // construct a small scanner and sinogram - shared_ptr scanner_sptr(new Scanner(Scanner::E953)); - scanner_sptr->set_num_rings(5); - shared_ptr proj_data_info_sptr(ProjDataInfo::ProjDataInfoCTI(scanner_sptr, + if (TOF_or_not) + { + std::cerr << "------ TOF data ----\n"; + shared_ptr scanner_sptr(new Scanner(Scanner::Discovery690)); + scanner_sptr->set_num_rings(4); + proj_data_info_sptr = std::move(ProjDataInfo::construct_proj_data_info(scanner_sptr, /*span=*/3, /*max_delta=*/4, /*num_views=*/16, - /*num_tang_poss=*/16)); + /*num_tang_poss=*/16, + /* arccorrected=*/false, + /* TOF_mash_factor=*/11)); + } + else + { + std::cerr << "------ non-TOF data ----\n"; + shared_ptr scanner_sptr(new Scanner(Scanner::E953)); + scanner_sptr->set_num_rings(5); + proj_data_info_sptr.reset(ProjDataInfo::ProjDataInfoCTI(scanner_sptr, + /*span=*/3, + /*max_delta=*/4, + /*num_views=*/16, + /*num_tang_poss=*/16)); + } shared_ptr exam_info_sptr(new ExamInfo(ImagingModality::PT)); proj_data_sptr.reset(new ProjDataInMemory(exam_info_sptr, proj_data_info_sptr)); for (int seg_num = proj_data_sptr->get_min_segment_num(); seg_num <= proj_data_sptr->get_max_segment_num(); ++seg_num) @@ -251,24 +270,23 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::construct_input_data // multiplicative term shared_ptr bin_norm_sptr(new TrivialBinNormalisation()); { + mult_proj_data_sptr.reset(new ProjDataInMemory(proj_data_sptr->get_exam_info_sptr(), - proj_data_sptr->get_proj_data_info_sptr()->create_shared_clone())); + proj_data_sptr->get_proj_data_info_sptr()->create_non_tof_clone())); for (int seg_num = proj_data_sptr->get_min_segment_num(); seg_num <= proj_data_sptr->get_max_segment_num(); ++seg_num) { - for (int timing_pos_num = proj_data_sptr->get_min_tof_pos_num(); timing_pos_num <= proj_data_sptr->get_max_tof_pos_num(); - ++timing_pos_num) - { - SegmentByView segment = proj_data_sptr->get_empty_segment_by_view(seg_num, false, timing_pos_num); - // fill in some crazy values - float value = 0; - for (SegmentByView::full_iterator iter = segment.begin_all(); iter != segment.end_all(); ++iter) - { - value = float(fabs(seg_num * value - .2)); // needs to be positive for Poisson - *iter = value; - } - segment /= 0.5F * segment.find_max(); // normalise to 2 (to avoid large numbers) - mult_proj_data_sptr->set_segment(segment); - } + { + auto segment = mult_proj_data_sptr->get_empty_segment_by_view(seg_num); + // fill in some crazy values + float value = 0; + for (auto iter = segment.begin_all(); iter != segment.end_all(); ++iter) + { + value = float(fabs(seg_num * value - .2)); // needs to be positive for Poisson + *iter = value; + } + segment /= 0.5F * segment.find_max(); // normalise to 2 (to avoid large numbers) + mult_proj_data_sptr->set_segment(segment); + } } bin_norm_sptr.reset(new BinNormalisationFromProjData(mult_proj_data_sptr)); } @@ -305,8 +323,8 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::construct_input_data shared_ptr proj_pair_sptr(new ProjectorByBinPairUsingProjMatrixByBin(proj_matrix_sptr)); objective_function.set_projector_pair_sptr(proj_pair_sptr); /* - void set_frame_num(const int); - void set_frame_definitions(const TimeFrameDefinitions&); + void set_frame_num(const int); + void set_frame_definitions(const TimeFrameDefinitions&); */ objective_function.set_normalisation_sptr(bin_norm_sptr); objective_function.set_additive_proj_data_sptr(add_proj_data_sptr); @@ -323,9 +341,18 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::run_tests() Verbosity::set(0); #if 1 - shared_ptr density_sptr; - construct_input_data(density_sptr); - this->run_tests_for_objective_function(*this->objective_function_sptr, *density_sptr); + { + shared_ptr density_sptr; + construct_input_data(density_sptr, /*TOF_or_not=*/false); + this->run_tests_for_objective_function(*this->objective_function_sptr, *density_sptr); + } + if (this->proj_data_filename == 0) + { + shared_ptr density_sptr; + construct_input_data(density_sptr, /*TOF_or_not=*/true); + this->run_tests_for_objective_function(*this->objective_function_sptr, *density_sptr); + } + #else // alternative that gets the objective function from an OSMAPOSL .par file // currently disabled