From e0f906f7b2f5749ee82fa09ea467d567f07cf57e Mon Sep 17 00:00:00 2001 From: NicoleJurjew Date: Thu, 14 Mar 2024 15:34:38 +0000 Subject: [PATCH 01/17] SSRB: add TOF mashing funcionality (BETA) --- src/buildblock/SSRB.cxx | 82 ++++++++++++++++++++++------------------- 1 file changed, 45 insertions(+), 37 deletions(-) diff --git a/src/buildblock/SSRB.cxx b/src/buildblock/SSRB.cxx index 7cd38d7ee8..0818b05813 100644 --- a/src/buildblock/SSRB.cxx +++ b/src/buildblock/SSRB.cxx @@ -47,8 +47,8 @@ SSRB(const ProjDataInfo& in_proj_data_info, const int max_in_segment_num_to_process_argument, const int num_tof_bins_to_combine) { - if (num_tof_bins_to_combine != 1) - error("SSRB: num_tof_bins_to_combine (%d) currently needs to be 1", num_tof_bins_to_combine); + // if (num_tof_bins_to_combine != 1) + // error("SSRB: num_tof_bins_to_combine (%d) currently needs to be 1", num_tof_bins_to_combine); if (num_segments_to_combine % 2 == 0) error("SSRB: num_segments_to_combine (%d) needs to be odd\n", num_segments_to_combine); const int max_in_segment_num_to_process = max_in_segment_num_to_process_argument >= 0 ? max_in_segment_num_to_process_argument @@ -165,8 +165,8 @@ SSRB(const string& output_filename, void SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) { - if (in_proj_data.get_proj_data_info_sptr()->is_tof_data()) - error("SSRB for TOF data is not currently implemented.\n"); + // if (in_proj_data.get_proj_data_info_sptr()->is_tof_data()) + // error("SSRB for TOF data is not currently implemented.\n"); const shared_ptr in_proj_data_info_sptr = dynamic_pointer_cast(in_proj_data.get_proj_data_info_sptr()); @@ -184,6 +184,7 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) } const int num_views_to_combine = in_proj_data.get_num_views() / out_proj_data.get_num_views(); + const int num_tof_bins_to_combine = in_proj_data.get_num_tof_poss() / out_proj_data.get_num_tof_poss(); if (in_proj_data.get_min_view_num() != 0 || out_proj_data.get_min_view_num() != 0) error("SSRB can only mash views when min_view_num==0\n"); @@ -238,47 +239,54 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) = out_proj_data.get_empty_sinogram(out_proj_data.get_min_axial_pos_num(out_segment_num), out_segment_num); Sinogram in_sino = in_proj_data.get_empty_sinogram(in_proj_data.get_min_axial_pos_num(out_segment_num), out_segment_num); - - for (int out_ax_pos_num = out_proj_data.get_min_axial_pos_num(out_segment_num); - out_ax_pos_num <= out_proj_data.get_max_axial_pos_num(out_segment_num); - ++out_ax_pos_num) + for (int out_timing_pos_num = out_proj_data.get_min_tof_pos_num(); + out_timing_pos_num <= out_proj_data.get_max_tof_pos_num(); + ++out_timing_pos_num) { - out_sino = out_proj_data.get_empty_sinogram(out_ax_pos_num, out_segment_num); - // get_m could be replaced by get_t - const float out_m = out_proj_data_info_sptr->get_m(Bin(out_segment_num, 0, out_ax_pos_num, 0)); + for (int out_ax_pos_num = out_proj_data.get_min_axial_pos_num(out_segment_num); + out_ax_pos_num <= out_proj_data.get_max_axial_pos_num(out_segment_num); + ++out_ax_pos_num) + { + out_sino = out_proj_data.get_empty_sinogram(out_ax_pos_num, out_segment_num, false, out_timing_pos_num); + // get_m could be replaced by get_t + const float out_m = out_proj_data_info_sptr->get_m(Bin(out_segment_num, 0, out_ax_pos_num, out_timing_pos_num)); - unsigned int num_in_ax_pos = 0; + unsigned int num_in_ax_pos = 0; - for (int in_segment_num = in_min_segment_num; in_segment_num <= in_max_segment_num; ++in_segment_num) - for (int in_ax_pos_num = in_proj_data.get_min_axial_pos_num(in_segment_num); - in_ax_pos_num <= in_proj_data.get_max_axial_pos_num(in_segment_num); - ++in_ax_pos_num) - { - const float in_m = in_proj_data_info_sptr->get_m(Bin(in_segment_num, 0, in_ax_pos_num, 0)); - if (fabs(out_m - in_m) < 1E-4) + for (int in_timing_pos_num = in_proj_data.get_min_tof_pos_num(); + in_timing_pos_num <= in_proj_data.get_max_tof_pos_num(); + ++in_timing_pos_num) + for (int in_segment_num = in_min_segment_num; in_segment_num <= in_max_segment_num; ++in_segment_num) + for (int in_ax_pos_num = in_proj_data.get_min_axial_pos_num(in_segment_num); + in_ax_pos_num <= in_proj_data.get_max_axial_pos_num(in_segment_num); + ++in_ax_pos_num) { - ++num_in_ax_pos; + const float in_m = in_proj_data_info_sptr->get_m(Bin(in_segment_num, 0, in_ax_pos_num, in_timing_pos_num)); + if (fabs(out_m - in_m) < 1E-4) + { + ++num_in_ax_pos; - in_sino = in_proj_data.get_sinogram(in_ax_pos_num, in_segment_num); - for (int in_view_num = in_proj_data.get_min_view_num(); in_view_num <= in_proj_data.get_max_view_num(); - ++in_view_num) - for (int tangential_pos_num - = max(in_proj_data.get_min_tangential_pos_num(), out_proj_data.get_min_tangential_pos_num()); - tangential_pos_num - <= min(in_proj_data.get_max_tangential_pos_num(), out_proj_data.get_max_tangential_pos_num()); - ++tangential_pos_num) - out_sino[in_view_num / num_views_to_combine][tangential_pos_num] - += in_sino[in_view_num][tangential_pos_num]; + in_sino = in_proj_data.get_sinogram(in_ax_pos_num, in_segment_num, false, in_timing_pos_num); + for (int in_view_num = in_proj_data.get_min_view_num(); in_view_num <= in_proj_data.get_max_view_num(); + ++in_view_num) + for (int tangential_pos_num + = max(in_proj_data.get_min_tangential_pos_num(), out_proj_data.get_min_tangential_pos_num()); + tangential_pos_num + <= min(in_proj_data.get_max_tangential_pos_num(), out_proj_data.get_max_tangential_pos_num()); + ++tangential_pos_num) + out_sino[in_view_num / num_views_to_combine][tangential_pos_num] + += in_sino[in_view_num][tangential_pos_num]; - break; // out of loop over ax_pos as we found where to put it + break; // out of loop over ax_pos as we found where to put it + } } - } - if (do_norm && num_in_ax_pos != 0) - out_sino /= static_cast(num_in_ax_pos * num_views_to_combine); - if (num_in_ax_pos == 0) - warning("SSRB: no sinograms contributing to output segment %d, ax_pos %d\n", out_segment_num, out_ax_pos_num); + if (do_norm && num_in_ax_pos != 0) + out_sino /= static_cast(num_in_ax_pos * num_views_to_combine); + if (num_in_ax_pos == 0) + warning("SSRB: no sinograms contributing to output segment %d, ax_pos %d\n", out_segment_num, out_ax_pos_num); - out_proj_data.set_sinogram(out_sino); + out_proj_data.set_sinogram(out_sino); + } } } } From 80d983b06f4f4fdb2c472b17bf1a35af226fcbac Mon Sep 17 00:00:00 2001 From: NicoleJurjew Date: Thu, 14 Mar 2024 15:37:26 +0000 Subject: [PATCH 02/17] SSRB: formatting (clang) --- src/buildblock/SSRB.cxx | 52 +++++++++++++++++++++-------------------- 1 file changed, 27 insertions(+), 25 deletions(-) diff --git a/src/buildblock/SSRB.cxx b/src/buildblock/SSRB.cxx index 0818b05813..8fadb572c3 100644 --- a/src/buildblock/SSRB.cxx +++ b/src/buildblock/SSRB.cxx @@ -253,33 +253,35 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) unsigned int num_in_ax_pos = 0; - for (int in_timing_pos_num = in_proj_data.get_min_tof_pos_num(); - in_timing_pos_num <= in_proj_data.get_max_tof_pos_num(); - ++in_timing_pos_num) - for (int in_segment_num = in_min_segment_num; in_segment_num <= in_max_segment_num; ++in_segment_num) - for (int in_ax_pos_num = in_proj_data.get_min_axial_pos_num(in_segment_num); - in_ax_pos_num <= in_proj_data.get_max_axial_pos_num(in_segment_num); - ++in_ax_pos_num) - { - const float in_m = in_proj_data_info_sptr->get_m(Bin(in_segment_num, 0, in_ax_pos_num, in_timing_pos_num)); - if (fabs(out_m - in_m) < 1E-4) - { - ++num_in_ax_pos; + for (int in_timing_pos_num = in_proj_data.get_min_tof_pos_num(); + in_timing_pos_num <= in_proj_data.get_max_tof_pos_num(); + ++in_timing_pos_num) + for (int in_segment_num = in_min_segment_num; in_segment_num <= in_max_segment_num; ++in_segment_num) + for (int in_ax_pos_num = in_proj_data.get_min_axial_pos_num(in_segment_num); + in_ax_pos_num <= in_proj_data.get_max_axial_pos_num(in_segment_num); + ++in_ax_pos_num) + { + const float in_m + = in_proj_data_info_sptr->get_m(Bin(in_segment_num, 0, in_ax_pos_num, in_timing_pos_num)); + if (fabs(out_m - in_m) < 1E-4) + { + ++num_in_ax_pos; - in_sino = in_proj_data.get_sinogram(in_ax_pos_num, in_segment_num, false, in_timing_pos_num); - for (int in_view_num = in_proj_data.get_min_view_num(); in_view_num <= in_proj_data.get_max_view_num(); - ++in_view_num) - for (int tangential_pos_num - = max(in_proj_data.get_min_tangential_pos_num(), out_proj_data.get_min_tangential_pos_num()); - tangential_pos_num - <= min(in_proj_data.get_max_tangential_pos_num(), out_proj_data.get_max_tangential_pos_num()); - ++tangential_pos_num) - out_sino[in_view_num / num_views_to_combine][tangential_pos_num] - += in_sino[in_view_num][tangential_pos_num]; + in_sino = in_proj_data.get_sinogram(in_ax_pos_num, in_segment_num, false, in_timing_pos_num); + for (int in_view_num = in_proj_data.get_min_view_num(); + in_view_num <= in_proj_data.get_max_view_num(); + ++in_view_num) + for (int tangential_pos_num + = max(in_proj_data.get_min_tangential_pos_num(), out_proj_data.get_min_tangential_pos_num()); + tangential_pos_num + <= min(in_proj_data.get_max_tangential_pos_num(), out_proj_data.get_max_tangential_pos_num()); + ++tangential_pos_num) + out_sino[in_view_num / num_views_to_combine][tangential_pos_num] + += in_sino[in_view_num][tangential_pos_num]; - break; // out of loop over ax_pos as we found where to put it - } - } + break; // out of loop over ax_pos as we found where to put it + } + } if (do_norm && num_in_ax_pos != 0) out_sino /= static_cast(num_in_ax_pos * num_views_to_combine); if (num_in_ax_pos == 0) From 854c5c68b2988bb8b0e05df2bd34adab2ed4312f Mon Sep 17 00:00:00 2001 From: NicoleJurjew Date: Fri, 15 Mar 2024 13:33:57 +0000 Subject: [PATCH 03/17] SSRB: include number of TOF bins to combine added in utility --- src/buildblock/SSRB.cxx | 4 ++-- src/include/stir/SSRB.h | 2 +- src/utilities/SSRB.cxx | 9 ++++++--- 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/buildblock/SSRB.cxx b/src/buildblock/SSRB.cxx index 8fadb572c3..1597ece862 100644 --- a/src/buildblock/SSRB.cxx +++ b/src/buildblock/SSRB.cxx @@ -47,8 +47,8 @@ SSRB(const ProjDataInfo& in_proj_data_info, const int max_in_segment_num_to_process_argument, const int num_tof_bins_to_combine) { - // if (num_tof_bins_to_combine != 1) - // error("SSRB: num_tof_bins_to_combine (%d) currently needs to be 1", num_tof_bins_to_combine); + if (in_proj_data_info.get_num_tof_poss() % num_tof_bins_to_combine == 0) + error("SSRB: num_tof_bins_to_combine needs to divide the number of TOF bins\n"); if (num_segments_to_combine % 2 == 0) error("SSRB: num_segments_to_combine (%d) needs to be odd\n", num_segments_to_combine); const int max_in_segment_num_to_process = max_in_segment_num_to_process_argument >= 0 ? max_in_segment_num_to_process_argument diff --git a/src/include/stir/SSRB.h b/src/include/stir/SSRB.h index 095bd6d8c1..926bc366d8 100644 --- a/src/include/stir/SSRB.h +++ b/src/include/stir/SSRB.h @@ -91,7 +91,7 @@ ProjDataInfo* SSRB(const ProjDataInfo& in_proj_data_info, Default value -1 means 'do all segments'. \param do_normalisation (default true) wether to normalise the output sinograms corresponding to how many input sinograms contribute to them. - \param num_tof_bins_to_combine currently has to be 1. If it doesn't, error() will be called. + \param num_tof_bins_to_combine defaults to 1, so TOF bins are not combined. \see SSRB(const ProjDataInfo& in_proj_data_info, const int num_segments_to_combine, diff --git a/src/utilities/SSRB.cxx b/src/utilities/SSRB.cxx index 8cde13ce22..0879f875a2 100644 --- a/src/utilities/SSRB.cxx +++ b/src/utilities/SSRB.cxx @@ -81,7 +81,7 @@ print_usage_and_exit(const std::string& prog_name) << prog_name << " [-t num_tangential_poss_to_trim] \\\n" << "\toutput_filename input_projdata_name \\\n" << "\t[num_segments_to_combine \\\n" - << "\t[ num_views_to_combine [do_norm [max_in_segment_num_to_process ]]]]\n" + << "\t[ num_views_to_combine [do_norm [max_in_segment_num_to_process [num_tof_bins_to_combine ]]]]]\n" << "num_segments_to_combine has to be odd. It is used as the number of segments\n" << " in the original data to combine.\n" << "num_views_to_combine has to be at least 1 (which is the default)\n" @@ -89,6 +89,7 @@ print_usage_and_exit(const std::string& prog_name) << " of tangential positions.\n" << "do_norm has to be 1 (normalise the result, which is the default) or 0\n" << "max_in_segment_num_to_process defaults to all segments\n" + << "num_tof_bins_to_combine defaults to 1, so TOF bins are not combined\n" << "\n\t - Template Based SSRB method\n" << prog_name << "--template template_projdata_filename output_filename input_projdata_name [do_norm]\n" << "template_projdata_filename the format of the output sinogram.\n" @@ -114,6 +115,7 @@ classic_SSRB(int argc, char** argv) const int num_views_to_combine = argc <= 4 ? 1 : atoi(argv[4]); const bool do_norm = argc <= 5 ? true : atoi(argv[5]) != 0; const int max_segment_num_to_process = argc <= 6 ? -1 : atoi(argv[6]); + const int num_tof_bins_to_combine = argc <= 7 ? 1 : atoi(argv[7]); // do standard SSRB SSRB(output_filename, *in_projdata_ptr, @@ -121,7 +123,8 @@ classic_SSRB(int argc, char** argv) num_views_to_combine, num_tangential_poss_to_trim, do_norm, - max_segment_num_to_process); + max_segment_num_to_process, + num_tof_bins_to_combine); } void @@ -140,7 +143,7 @@ template_based_SSRB(int argc, char** argv) int main(int argc, char** argv) { - if (argc > 7 || argc < 3) + if (argc > 8 || argc < 3) { print_usage_and_exit(argv[0]); } From b0001ab812e747c8e0069a4bebce25e38df41eff Mon Sep 17 00:00:00 2001 From: NicoleJurjew Date: Fri, 15 Mar 2024 15:21:34 +0000 Subject: [PATCH 04/17] remove check for tof-bins --- src/buildblock/SSRB.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/buildblock/SSRB.cxx b/src/buildblock/SSRB.cxx index 1597ece862..f3863baac0 100644 --- a/src/buildblock/SSRB.cxx +++ b/src/buildblock/SSRB.cxx @@ -47,8 +47,8 @@ SSRB(const ProjDataInfo& in_proj_data_info, const int max_in_segment_num_to_process_argument, const int num_tof_bins_to_combine) { - if (in_proj_data_info.get_num_tof_poss() % num_tof_bins_to_combine == 0) - error("SSRB: num_tof_bins_to_combine needs to divide the number of TOF bins\n"); + // if (in_proj_data_info.get_num_tof_poss() % num_tof_bins_to_combine == 0) + // error("SSRB: num_tof_bins_to_combine needs to divide the number of TOF bins\n"); if (num_segments_to_combine % 2 == 0) error("SSRB: num_segments_to_combine (%d) needs to be odd\n", num_segments_to_combine); const int max_in_segment_num_to_process = max_in_segment_num_to_process_argument >= 0 ? max_in_segment_num_to_process_argument From 53a0446b1d09a84fa2a16c5fa8fd16c0125debd5 Mon Sep 17 00:00:00 2001 From: NicoleJurjew Date: Mon, 3 Jun 2024 11:45:13 +0100 Subject: [PATCH 05/17] src/buildblock/SSRB.cxx: bugfix: connect in_tof_bin to out_tof_bin src/utilities/SSRB.cxx: formatting of help doc-string --- src/buildblock/SSRB.cxx | 13 +++++++++++++ src/utilities/SSRB.cxx | 31 ++++++++++++++++--------------- 2 files changed, 29 insertions(+), 15 deletions(-) diff --git a/src/buildblock/SSRB.cxx b/src/buildblock/SSRB.cxx index f3863baac0..1f8642f072 100644 --- a/src/buildblock/SSRB.cxx +++ b/src/buildblock/SSRB.cxx @@ -253,6 +253,12 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) unsigned int num_in_ax_pos = 0; + Bin out_bin(out_segment_num, 0, out_ax_pos_num, 0, out_timing_pos_num); + + // get edges of TOF bin, currently only exposed via sampling + const float out_lower_k = out_proj_data_info_sptr->get_k(out_bin) - out_proj_data_info_sptr->get_sampling_in_k(out_bin)/2; + const float out_higher_k = out_proj_data_info_sptr->get_k(out_bin) + out_proj_data_info_sptr->get_sampling_in_k(out_bin)/2; + for (int in_timing_pos_num = in_proj_data.get_min_tof_pos_num(); in_timing_pos_num <= in_proj_data.get_max_tof_pos_num(); ++in_timing_pos_num) @@ -268,6 +274,13 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) ++num_in_ax_pos; in_sino = in_proj_data.get_sinogram(in_ax_pos_num, in_segment_num, false, in_timing_pos_num); + + Bin in_bin(in_segment_num, 0, in_ax_pos_num, 0, in_timing_pos_num); + const float in_k = in_proj_data_info_sptr->get_k(in_bin); + + // check if in_timing_pos_num is in the range for the out bin or not + if (in_k < out_lower_k || in_k >= out_higher_k) + continue; for (int in_view_num = in_proj_data.get_min_view_num(); in_view_num <= in_proj_data.get_max_view_num(); ++in_view_num) diff --git a/src/utilities/SSRB.cxx b/src/utilities/SSRB.cxx index 0879f875a2..8b51dae537 100644 --- a/src/utilities/SSRB.cxx +++ b/src/utilities/SSRB.cxx @@ -77,24 +77,25 @@ print_usage_and_exit(const std::string& prog_name) { cerr << "Usage:\n\n" << "Two options:\n" - << "\t - Classical STIR SSRB method\n" + << "Classical STIR SSRB method:\n" << prog_name << " [-t num_tangential_poss_to_trim] \\\n" << "\toutput_filename input_projdata_name \\\n" << "\t[num_segments_to_combine \\\n" - << "\t[ num_views_to_combine [do_norm [max_in_segment_num_to_process [num_tof_bins_to_combine ]]]]]\n" - << "num_segments_to_combine has to be odd. It is used as the number of segments\n" - << " in the original data to combine.\n" - << "num_views_to_combine has to be at least 1 (which is the default)\n" - << "num_tangential_poss_to_trim has to be smaller than the available number\n" - << " of tangential positions.\n" - << "do_norm has to be 1 (normalise the result, which is the default) or 0\n" - << "max_in_segment_num_to_process defaults to all segments\n" - << "num_tof_bins_to_combine defaults to 1, so TOF bins are not combined\n" - << "\n\t - Template Based SSRB method\n" - << prog_name << "--template template_projdata_filename output_filename input_projdata_name [do_norm]\n" - << "template_projdata_filename the format of the output sinogram.\n" - << "output_filename, input_projdata_name, and do_norm are as above.\n" - << "SSRB act in the same way as Classical but allows for even num_segments_to_combine.\n"; + << "\t[num_views_to_combine [do_norm [max_in_segment_num_to_process\\\n" + << "\t[num_tof_bins_to_combine ]]]]]\n" + << "\t-num_segments_to_combine has to be odd. It is used as the number of segments\n" + << "\t in the original data to combine.\n" + << "\t-num_views_to_combine has to be at least 1 (which is the default)\n" + << "\t-num_tangential_poss_to_trim has to be smaller than the available number\n" + << "\t of tangential positions.\n" + << "\t-do_norm has to be 1 (normalise the result, which is the default) or 0\n" + << "\t-max_in_segment_num_to_process defaults to all segments\n" + << "\t-num_tof_bins_to_combine defaults to 1, so TOF bins are not combined\n" + << "\nTemplate Based SSRB method:\n" + << prog_name << " --template template_projdata_filename output_filename input_projdata_name [do_norm]\n" + << "\t-template_projdata_filename the format of the output sinogram.\n" + << "\t-output_filename, input_projdata_name, and do_norm are as above.\n" + << "\tSSRB acts in the same way as Classical but allows for even num_segments_to_combine.\n"; exit(EXIT_FAILURE); } From 884a65ef61e42eff34565ed81b23883ef2c69d40 Mon Sep 17 00:00:00 2001 From: NicoleJurjew Date: Mon, 3 Jun 2024 11:50:17 +0100 Subject: [PATCH 06/17] src/buildblock/SSRB.cxx: apply clang format --- src/buildblock/SSRB.cxx | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/buildblock/SSRB.cxx b/src/buildblock/SSRB.cxx index 1f8642f072..36423de835 100644 --- a/src/buildblock/SSRB.cxx +++ b/src/buildblock/SSRB.cxx @@ -256,8 +256,10 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) Bin out_bin(out_segment_num, 0, out_ax_pos_num, 0, out_timing_pos_num); // get edges of TOF bin, currently only exposed via sampling - const float out_lower_k = out_proj_data_info_sptr->get_k(out_bin) - out_proj_data_info_sptr->get_sampling_in_k(out_bin)/2; - const float out_higher_k = out_proj_data_info_sptr->get_k(out_bin) + out_proj_data_info_sptr->get_sampling_in_k(out_bin)/2; + const float out_lower_k + = out_proj_data_info_sptr->get_k(out_bin) - out_proj_data_info_sptr->get_sampling_in_k(out_bin) / 2; + const float out_higher_k + = out_proj_data_info_sptr->get_k(out_bin) + out_proj_data_info_sptr->get_sampling_in_k(out_bin) / 2; for (int in_timing_pos_num = in_proj_data.get_min_tof_pos_num(); in_timing_pos_num <= in_proj_data.get_max_tof_pos_num(); @@ -274,13 +276,13 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) ++num_in_ax_pos; in_sino = in_proj_data.get_sinogram(in_ax_pos_num, in_segment_num, false, in_timing_pos_num); - + Bin in_bin(in_segment_num, 0, in_ax_pos_num, 0, in_timing_pos_num); const float in_k = in_proj_data_info_sptr->get_k(in_bin); - + // check if in_timing_pos_num is in the range for the out bin or not if (in_k < out_lower_k || in_k >= out_higher_k) - continue; + continue; for (int in_view_num = in_proj_data.get_min_view_num(); in_view_num <= in_proj_data.get_max_view_num(); ++in_view_num) From 73f1139552bca7faa8d1c00537c7f36c624795a7 Mon Sep 17 00:00:00 2001 From: NicoleJurjew Date: Mon, 3 Jun 2024 12:00:28 +0100 Subject: [PATCH 07/17] src/swig/stir.i & src/swig/stir_projdata.i: add iSSRB and SSRB to the SWIG-interface --- src/swig/stir.i | 3 +++ src/swig/stir_projdata.i | 12 ++++++++++++ 2 files changed, 15 insertions(+) diff --git a/src/swig/stir.i b/src/swig/stir.i index 9ac439079e..3966306cba 100644 --- a/src/swig/stir.i +++ b/src/swig/stir.i @@ -162,6 +162,9 @@ #include "stir/scatter/SingleScatterSimulation.h" #include "stir/scatter/CreateTailMaskFromACFs.h" +#include "stir/SSRB.h" +#include "stir/inverse_SSRB.h" + #include #include #include diff --git a/src/swig/stir_projdata.i b/src/swig/stir_projdata.i index 6a1e0eb433..25d534e683 100644 --- a/src/swig/stir_projdata.i +++ b/src/swig/stir_projdata.i @@ -64,6 +64,7 @@ ADD_REPR(stir::DetectionPositionPair, %arg(*$self)) %include "stir/ViewgramIndices.h" %include "stir/SinogramIndices.h" %include "stir/Bin.h" + ADD_REPR(stir::Bin, %arg(*$self)) @@ -274,3 +275,14 @@ namespace stir { //%template(SharedProjData) boost::shared_ptr; } + + +%include "stir/SSRB.h" +%include "stir/inverse_SSRB.h" + +%newobject stir::SSRB; /*(const stir::ProjDataInfo& in_proj_data_info, + const int num_segments_to_combine, + const int num_views_to_combine = 1, + const int num_tangential_poss_to_trim = 0, + const int max_in_segment_num_to_process = -1, + const int num_tof_bins_to_combine = 1);*/ \ No newline at end of file From 8948c7a46b93d60be9902bc2cab0eafed96f0c46 Mon Sep 17 00:00:00 2001 From: NicoleJurjew Date: Thu, 20 Jun 2024 11:51:41 +0100 Subject: [PATCH 08/17] src/buildblock/SSRB.cxx: corrected normalization --- src/buildblock/SSRB.cxx | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/buildblock/SSRB.cxx b/src/buildblock/SSRB.cxx index 36423de835..8256d19400 100644 --- a/src/buildblock/SSRB.cxx +++ b/src/buildblock/SSRB.cxx @@ -251,7 +251,7 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) // get_m could be replaced by get_t const float out_m = out_proj_data_info_sptr->get_m(Bin(out_segment_num, 0, out_ax_pos_num, out_timing_pos_num)); - unsigned int num_in_ax_pos = 0; + unsigned int num_in_sinos = 0; Bin out_bin(out_segment_num, 0, out_ax_pos_num, 0, out_timing_pos_num); @@ -273,16 +273,14 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) = in_proj_data_info_sptr->get_m(Bin(in_segment_num, 0, in_ax_pos_num, in_timing_pos_num)); if (fabs(out_m - in_m) < 1E-4) { - ++num_in_ax_pos; - - in_sino = in_proj_data.get_sinogram(in_ax_pos_num, in_segment_num, false, in_timing_pos_num); - Bin in_bin(in_segment_num, 0, in_ax_pos_num, 0, in_timing_pos_num); const float in_k = in_proj_data_info_sptr->get_k(in_bin); // check if in_timing_pos_num is in the range for the out bin or not if (in_k < out_lower_k || in_k >= out_higher_k) continue; + in_sino = in_proj_data.get_sinogram(in_ax_pos_num, in_segment_num, false, in_timing_pos_num); + ++num_in_sinos; for (int in_view_num = in_proj_data.get_min_view_num(); in_view_num <= in_proj_data.get_max_view_num(); ++in_view_num) @@ -297,9 +295,9 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) break; // out of loop over ax_pos as we found where to put it } } - if (do_norm && num_in_ax_pos != 0) - out_sino /= static_cast(num_in_ax_pos * num_views_to_combine); - if (num_in_ax_pos == 0) + if (do_norm && num_in_sinos != 0) + out_sino /= static_cast(num_in_sinos * num_views_to_combine); + if (num_in_sinos == 0) warning("SSRB: no sinograms contributing to output segment %d, ax_pos %d\n", out_segment_num, out_ax_pos_num); out_proj_data.set_sinogram(out_sino); From cc787f514efc33d2bbe860c4ae9561d8d2ba49b4 Mon Sep 17 00:00:00 2001 From: NicoleJurjew Date: Mon, 24 Jun 2024 10:22:51 +0100 Subject: [PATCH 09/17] add TOF mashing to release notes --- documentation/release_6.2.htm | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/documentation/release_6.2.htm b/documentation/release_6.2.htm index ccab3e59ea..14cb81786a 100644 --- a/documentation/release_6.2.htm +++ b/documentation/release_6.2.htm @@ -78,6 +78,11 @@

Changed functionality


PR #1458 +
  • + SSRB now allows to mash TOF bins. +
    + PR #1464 +
    • From 2493ae0e8243bf322458f6685c08643f9e38268a Mon Sep 17 00:00:00 2001 From: NicoleJurjew Date: Mon, 24 Jun 2024 14:08:43 +0100 Subject: [PATCH 10/17] fix SSRB for nonTOF data --- src/buildblock/SSRB.cxx | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/buildblock/SSRB.cxx b/src/buildblock/SSRB.cxx index 8256d19400..93907eed0e 100644 --- a/src/buildblock/SSRB.cxx +++ b/src/buildblock/SSRB.cxx @@ -256,10 +256,16 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) Bin out_bin(out_segment_num, 0, out_ax_pos_num, 0, out_timing_pos_num); // get edges of TOF bin, currently only exposed via sampling + // for non-TOF data, the sampling in k is 0, therefore condition below is never met + // therefore: for non-TOF set out_lower_k to -1E20F and out_higher_k to 1E20F const float out_lower_k - = out_proj_data_info_sptr->get_k(out_bin) - out_proj_data_info_sptr->get_sampling_in_k(out_bin) / 2; + = out_proj_data_info_sptr->is_tof_data() + ? (out_proj_data_info_sptr->get_k(out_bin) - out_proj_data_info_sptr->get_sampling_in_k(out_bin) / 2) + : -1E20F; const float out_higher_k - = out_proj_data_info_sptr->get_k(out_bin) + out_proj_data_info_sptr->get_sampling_in_k(out_bin) / 2; + = out_proj_data_info_sptr->is_tof_data() + ? (out_proj_data_info_sptr->get_k(out_bin) + out_proj_data_info_sptr->get_sampling_in_k(out_bin) / 2) + : 1E20F; for (int in_timing_pos_num = in_proj_data.get_min_tof_pos_num(); in_timing_pos_num <= in_proj_data.get_max_tof_pos_num(); From ba80e3b23542ba337fcf620a5f8877dec8d289e2 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Mon, 1 Jul 2024 16:05:01 +0100 Subject: [PATCH 11/17] use in/out_bin consistently This results in shorter code, but also consistent. There was a "problem" in the previous version that the TOF bin was actually used as tang_pos_num when determining get_m() (but it didn't create a problem as currently get_m() in independent of both). --- src/buildblock/SSRB.cxx | 27 ++++++++++----------------- 1 file changed, 10 insertions(+), 17 deletions(-) diff --git a/src/buildblock/SSRB.cxx b/src/buildblock/SSRB.cxx index 93907eed0e..3c1f4dd2b5 100644 --- a/src/buildblock/SSRB.cxx +++ b/src/buildblock/SSRB.cxx @@ -2,7 +2,7 @@ // /* Copyright (C) 2002- 2013, Hammersmith Imanet Ltd - Copyright (C) 2021, University College London + Copyright (C) 2021, 2024, University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 @@ -16,7 +16,7 @@ \brief Implementation of stir::SSRB \author Kris Thielemans - + \author Nicole Jurjew */ #include "stir/ProjDataFromStream.h" #include "stir/ProjDataInterfile.h" @@ -165,9 +165,6 @@ SSRB(const string& output_filename, void SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) { - // if (in_proj_data.get_proj_data_info_sptr()->is_tof_data()) - // error("SSRB for TOF data is not currently implemented.\n"); - const shared_ptr in_proj_data_info_sptr = dynamic_pointer_cast(in_proj_data.get_proj_data_info_sptr()); if (is_null_ptr(in_proj_data_info_sptr)) @@ -184,7 +181,6 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) } const int num_views_to_combine = in_proj_data.get_num_views() / out_proj_data.get_num_views(); - const int num_tof_bins_to_combine = in_proj_data.get_num_tof_poss() / out_proj_data.get_num_tof_poss(); if (in_proj_data.get_min_view_num() != 0 || out_proj_data.get_min_view_num() != 0) error("SSRB can only mash views when min_view_num==0\n"); @@ -247,17 +243,15 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) out_ax_pos_num <= out_proj_data.get_max_axial_pos_num(out_segment_num); ++out_ax_pos_num) { - out_sino = out_proj_data.get_empty_sinogram(out_ax_pos_num, out_segment_num, false, out_timing_pos_num); - // get_m could be replaced by get_t - const float out_m = out_proj_data_info_sptr->get_m(Bin(out_segment_num, 0, out_ax_pos_num, out_timing_pos_num)); + Bin out_bin(out_segment_num, 0, out_ax_pos_num, 0, out_timing_pos_num); + out_sino = out_proj_data.get_empty_sinogram(out_bin); + const float out_m = out_proj_data_info_sptr->get_m(out_bin); unsigned int num_in_sinos = 0; - Bin out_bin(out_segment_num, 0, out_ax_pos_num, 0, out_timing_pos_num); - // get edges of TOF bin, currently only exposed via sampling - // for non-TOF data, the sampling in k is 0, therefore condition below is never met - // therefore: for non-TOF set out_lower_k to -1E20F and out_higher_k to 1E20F + // for non-TOF data, the sampling in k is 0, which is incorrect and would lead to the TOF condition below never + // being met. Therefore: for non-TOF set out_lower_k to -1E20F and out_higher_k to 1E20F const float out_lower_k = out_proj_data_info_sptr->is_tof_data() ? (out_proj_data_info_sptr->get_k(out_bin) - out_proj_data_info_sptr->get_sampling_in_k(out_bin) / 2) @@ -275,17 +269,16 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) in_ax_pos_num <= in_proj_data.get_max_axial_pos_num(in_segment_num); ++in_ax_pos_num) { - const float in_m - = in_proj_data_info_sptr->get_m(Bin(in_segment_num, 0, in_ax_pos_num, in_timing_pos_num)); + Bin in_bin(in_segment_num, 0, in_ax_pos_num, 0, in_timing_pos_num); + const float in_m = in_proj_data_info_sptr->get_m(in_bin); if (fabs(out_m - in_m) < 1E-4) { - Bin in_bin(in_segment_num, 0, in_ax_pos_num, 0, in_timing_pos_num); const float in_k = in_proj_data_info_sptr->get_k(in_bin); // check if in_timing_pos_num is in the range for the out bin or not if (in_k < out_lower_k || in_k >= out_higher_k) continue; - in_sino = in_proj_data.get_sinogram(in_ax_pos_num, in_segment_num, false, in_timing_pos_num); + in_sino = in_proj_data.get_sinogram(in_bin); ++num_in_sinos; for (int in_view_num = in_proj_data.get_min_view_num(); in_view_num <= in_proj_data.get_max_view_num(); From d874fbfd3ee7a04dc5cd13d36f398b541381d8a4 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 3 Jul 2024 00:48:52 +0100 Subject: [PATCH 12/17] SSRB: norm should not count TOF bins even if do_norm=1, we should not normalise with the number of TOF bins. This is because the projectors will take the width of the TOF bin properly into account. (In contrast, for "span", they don't). --- src/buildblock/SSRB.cxx | 17 ++++++++++------- src/include/stir/SSRB.h | 37 +++++++++++++++++++++++-------------- 2 files changed, 33 insertions(+), 21 deletions(-) diff --git a/src/buildblock/SSRB.cxx b/src/buildblock/SSRB.cxx index 3c1f4dd2b5..2a2c46df65 100644 --- a/src/buildblock/SSRB.cxx +++ b/src/buildblock/SSRB.cxx @@ -247,7 +247,7 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) out_sino = out_proj_data.get_empty_sinogram(out_bin); const float out_m = out_proj_data_info_sptr->get_m(out_bin); - unsigned int num_in_sinos = 0; + unsigned int num_in_ax_pos = 0; // get edges of TOF bin, currently only exposed via sampling // for non-TOF data, the sampling in k is 0, which is incorrect and would lead to the TOF condition below never @@ -273,13 +273,15 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) const float in_m = in_proj_data_info_sptr->get_m(in_bin); if (fabs(out_m - in_m) < 1E-4) { - const float in_k = in_proj_data_info_sptr->get_k(in_bin); + if (in_timing_pos_num == in_proj_data.get_min_tof_pos_num()) + ++num_in_ax_pos; // count number of in_sinograms contributing (ignoring TOF) + const float in_k = in_proj_data_info_sptr->get_k(in_bin); // check if in_timing_pos_num is in the range for the out bin or not if (in_k < out_lower_k || in_k >= out_higher_k) continue; in_sino = in_proj_data.get_sinogram(in_bin); - ++num_in_sinos; + for (int in_view_num = in_proj_data.get_min_view_num(); in_view_num <= in_proj_data.get_max_view_num(); ++in_view_num) @@ -294,10 +296,11 @@ SSRB(ProjData& out_proj_data, const ProjData& in_proj_data, const bool do_norm) break; // out of loop over ax_pos as we found where to put it } } - if (do_norm && num_in_sinos != 0) - out_sino /= static_cast(num_in_sinos * num_views_to_combine); - if (num_in_sinos == 0) - warning("SSRB: no sinograms contributing to output segment %d, ax_pos %d\n", out_segment_num, out_ax_pos_num); + if (do_norm && num_in_ax_pos != 0) + out_sino /= static_cast(num_in_ax_pos * num_views_to_combine); + if (num_in_ax_pos == 0) + warning("SSRB: no sinograms contributing to output segment " + std::to_string(out_segment_num) + ", ax_pos " + + std::to_string(out_ax_pos_num) + ", tof_pos_num " + std::to_string(out_timing_pos_num)); out_proj_data.set_sinogram(out_sino); } diff --git a/src/include/stir/SSRB.h b/src/include/stir/SSRB.h index 926bc366d8..b2c13c80d8 100644 --- a/src/include/stir/SSRB.h +++ b/src/include/stir/SSRB.h @@ -11,7 +11,7 @@ */ /* Copyright (C) 2002- 2009, Hammersmith Imanet Ltd - Copyright (C) 2021, University College London + Copyright (C) 2021, 2024 University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 @@ -55,15 +55,15 @@ class ProjDataInfo; segments. This essentially increases the axial compression (or span in CTI terminology), see the STIR Glossary on axial compression. In addition, SSRB can introduce extra mashing (see the STIR Glossary) of the data, i.e. add - views together. + views together. Finally, it can also be used to combine TOF bins together. Here is how to determine which bins are discarded when trimming is used. For a certain num_tangential_poss, the range is from \verbatim -(num_tangential_poss/2) to -(num_tangential_poss/2) + num_tangential_poss - 1. \endverbatim - The new num_tangential_poss is simply set to old_num_tangential_poss - - num_tang_poss_to_trim. Note that because of this, if num_tang_poss_to_trim is + The new \c num_tangential_poss is simply set to \c old_num_tangential_poss - + \a num_tang_poss_to_trim. Note that because of this, if \a num_tang_poss_to_trim is negative, more (zero) bins will be added. \warning in_proj_data_info has to be (at least) of type ProjDataInfoCylindrical @@ -93,13 +93,9 @@ ProjDataInfo* SSRB(const ProjDataInfo& in_proj_data_info, corresponding to how many input sinograms contribute to them. \param num_tof_bins_to_combine defaults to 1, so TOF bins are not combined. - \see SSRB(const ProjDataInfo& in_proj_data_info, - const int num_segments_to_combine, - const int num_views_to_combine, - const int num_tang_poss_to_trim, - const int max_in_segment_num_to_process, - const int num_tof_bins_to_combine - ) for restrictions + \see SSRB(ProjData& out_projdata, + const ProjData& in_projdata, + const bool do_normalisation = true) */ void SSRB(const std::string& output_filename, const ProjData& in_projdata, @@ -118,10 +114,23 @@ void SSRB(const std::string& output_filename, ProjData::set_sinogram(). \param in_projdata input data \param do_normalisation (default true) wether to normalise the output sinograms - corresponding to how many input sinograms contribute to them. + corresponding to how many (ignoring TOF) input sinograms contribute to them. + Note that we do not normalise according to the number of TOF bins. + This is because the projectors will take the width of the TOF bin + properly into account (by integration over the TOF kernel). (In contrast, in spatial + direction, projectors are outputting "normalised" data, i.e. corresponding to the + line integral). - \warning in_proj_data_info has to be (at least) of type ProjDataInfoCylindrical - \warning TOF info has to match currently. If it doesn't, error() will be called. + + \warning \a in_projdata has to be (at least) of type ProjDataInfoCylindrical + + \see SSRB(const ProjDataInfo& in_proj_data_info, + const int num_segments_to_combine, + const int num_views_to_combine, + const int num_tang_poss_to_trim, + const int max_in_segment_num_to_process, + const int num_tof_bins_to_combine + ) for information on the rebinning. */ void SSRB(ProjData& out_projdata, const ProjData& in_projdata, const bool do_normalisation = true); From 34c40c8d06016abb841307c8a5dec869810ff95e Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 3 Jul 2024 08:30:27 +0100 Subject: [PATCH 13/17] make simulate_PET_data_for_tests.sh more flexible add a few arguments and configurable environment variables to be able to modify the output of the simulation --- .../simulate_PET_data_for_tests.sh | 61 ++++++++++++------- 1 file changed, 40 insertions(+), 21 deletions(-) diff --git a/recon_test_pack/simulate_PET_data_for_tests.sh b/recon_test_pack/simulate_PET_data_for_tests.sh index e261877eb6..94e884a0bf 100755 --- a/recon_test_pack/simulate_PET_data_for_tests.sh +++ b/recon_test_pack/simulate_PET_data_for_tests.sh @@ -1,13 +1,13 @@ -#! /bin/sh +#! /bin/sh -vx # A script to simulate some data used by other tests. # Careful: run_scatter_tests.sh compares with previously generated data -# so you cannot simply modify this one with adjusting that as well. +# so you cannot simply modify this one without adjusting that as well. # # This script is not intended to be run on its own! # # Copyright (C) 2011 - 2011-01-14, Hammersmith Imanet Ltd # Copyright (C) 2011-07-01 - 2011, Kris Thielemans -# Copyright (C) 2014, 2020, 2022 University College London +# Copyright (C) 2014, 2020, 2022, 2024 University College London # This file is part of STIR. # # SPDX-License-Identifier: Apache-2.0 @@ -25,6 +25,7 @@ echo command -v generate_image >/dev/null 2>&1 || { echo "generate_image not found or not executable. Aborting." >&2; exit 1; } echo "Using `command -v generate_image`" +generate_images=1 force_zero_view_offset=0 TOF=0 suffix="" @@ -46,10 +47,19 @@ do then suffix="$2" shift 1 + elif test "$1" = "--keep_images" + then + generate_images=0 elif test "$1" = "--help" then - echo "Usage: `basename $0` [--force_zero_view_offset] [--suffix sometext] [install_dir]" + echo "Usage: `basename $0` [--force_zero_view_offset] [--suffix sometext] [--keep_images] [install_dir]" echo "(where [] means that an argument is optional)" + echo "The suffix is appended to filenames (before adding the extensions .hs/.s)." + echo "--keep-images can be used to avoid regenerating the images (if you know they are correct)." + echo "You can choose data sizes different from the defaults by setting the following environment variables:" + echo " view_mash, span, max_rd, tof_mash" + echo "Randoms+scatter are set to a constant proj_data. you can set that value via the environment variable:" + echo " background_value" exit 1 else echo Warning: Unknown option "$1" @@ -71,38 +81,46 @@ fi LC_ALL=C export LC_ALL +if [ "$generate_images" -eq 1 ]; then + echo "=== make emission image" + generate_image generate_uniform_cylinder.par + echo "=== make attenuation image" + generate_image generate_atten_cylinder.par +fi -echo "=== make emission image" -generate_image generate_uniform_cylinder.par -echo "=== make attenuation image" -generate_image generate_atten_cylinder.par if [ "$TOF" -eq 0 ]; then - echo "=== create template sinogram (DSTE in 3D with max ring diff 2 to save time)" - template_sino=my_DSTE_3D_rd3_template.hs + : ${view_mash:=1} + : ${span:=2} + : ${max_rd:=2} + echo "=== create template sinogram (DSTE with view_mash=${view_mash}, max_ring_diff=${max_rd})" + template_sino=my_DSTE_3D_vm${view_mash}_span${span}_rd${max_rd}_template.hs cat > my_input.txt < my_input.txt < my_create_${template_sino}.log 2>&1 if [ $? -ne 0 ]; then echo "ERROR running create_projdata_template. Check my_create_${template_sino}.log"; exit 1; @@ -125,7 +143,8 @@ if [ $force_zero_view_offset -eq 1 ]; then fi # create sinograms -./simulate_data.sh my_uniform_cylinder.hv my_atten_image.hv ${template_sino} 10 ${suffix} +: ${background_value:=10} +./simulate_data.sh my_uniform_cylinder.hv my_atten_image.hv ${template_sino} ${background_value} ${suffix} if [ $? -ne 0 ]; then echo "Error running simulation" exit 1 From d96d8fcf03b672c5b0cff0489691699eba28d5b7 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 3 Jul 2024 14:38:14 +0100 Subject: [PATCH 14/17] avoid caching matrix for forward projection Using the cache has very little impact on run-time for a single forward projection, but affects memory a lot. This is especially important for TOF. --- recon_test_pack/forward_projector_proj_matrix_ray_tracing.par | 1 + 1 file changed, 1 insertion(+) diff --git a/recon_test_pack/forward_projector_proj_matrix_ray_tracing.par b/recon_test_pack/forward_projector_proj_matrix_ray_tracing.par index 84c7c2d5d6..415a33f42f 100644 --- a/recon_test_pack/forward_projector_proj_matrix_ray_tracing.par +++ b/recon_test_pack/forward_projector_proj_matrix_ray_tracing.par @@ -4,6 +4,7 @@ type:=Matrix Forward Projector Using Matrix Parameters := Matrix type := Ray Tracing Ray tracing matrix parameters := + disable caching:= 1 number of rays in tangential direction to trace for each bin := 10 End Ray tracing matrix parameters := End Forward Projector Using Matrix Parameters := From f9cbda344c66d6be347ace291ba8c5cdfc033cc6 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 3 Jul 2024 14:51:38 +0100 Subject: [PATCH 15/17] add recon_test_pack/run_test_SSRB.sh --- documentation/release_6.2.htm | 11 ++ recon_test_pack/README.txt | 5 +- recon_test_pack/run_test_SSRB.sh | 197 +++++++++++++++++++++++++++++++ 3 files changed, 212 insertions(+), 1 deletion(-) create mode 100755 recon_test_pack/run_test_SSRB.sh diff --git a/documentation/release_6.2.htm b/documentation/release_6.2.htm index 14cb81786a..f7600fa0d3 100644 --- a/documentation/release_6.2.htm +++ b/documentation/release_6.2.htm @@ -169,6 +169,17 @@

    C++ tests

    +

    recon_test_pack

    +
      +
    • + The output of simulate_PET_data_for_tests.sh can now varied by setting environment variables, e.g. + max_rd. +
      • +
    • + New test run_test_SSRB.sh +
      • +
    + diff --git a/recon_test_pack/README.txt b/recon_test_pack/README.txt index 0a854c4a9c..8572da145f 100644 --- a/recon_test_pack/README.txt +++ b/recon_test_pack/README.txt @@ -1,5 +1,5 @@ -README file for STIR recon_test_pack version 6.1 (and later versions?) +README file for STIR recon_test_pack version 6.2 (and later versions?) ---------------------------------------------------------------------- This test pack runs some simple tests to check if various STIR reconstruction @@ -85,6 +85,9 @@ Similarly, a test for motion corrected reconstruction of gated data is run as sh run_test_simulate_and_recon_with_motion.sh [ --mpicmd cmd] [optional_install_path] +Finally, there are some tests on using SSRB to rebin the data to lower dimensions + +sh run_test_SSRB.sh [ --mpicmd cmd] [optional_install_path] Testing SPECT reconstructions diff --git a/recon_test_pack/run_test_SSRB.sh b/recon_test_pack/run_test_SSRB.sh new file mode 100755 index 0000000000..627ad04025 --- /dev/null +++ b/recon_test_pack/run_test_SSRB.sh @@ -0,0 +1,197 @@ +#! /bin/sh +# A script to check to see if reconstruction of simulated data gives the expected result. +# +# Copyright (C) 2011 - 2011-01-14, Hammersmith Imanet Ltd +# Copyright (C) 2011-07-01 - 2011, Kris Thielemans +# Copyright (C) 2014, 2022, 2024 University College London +# This file is part of STIR. +# +# SPDX-License-Identifier: Apache-2.0 +# +# See STIR/LICENSE.txt for details +# +# Author Kris Thielemans +# + +echo This script should work with STIR version 6.2. If you have +echo a later version, you might have to update your test pack. +echo Please check the web site. +echo + +# +# Options +# +MPIRUN="" +# +# Parse option arguments (--) +# Note that the -- is required to suppress interpretation of $1 as options +# to expr +# +while test `expr -- "$1" : "--.*"` -gt 0 +do + + if test "$1" = "--mpicmd" + then + MPIRUN="$2" + shift 1 + elif test "$1" = "--help" + then + echo "Usage: `basename $0` [--mpicmd somecmd] [install_dir]" + echo "(where [] means that an argument is optional)" + echo "See README.txt for more info." + exit 1 + else + echo Warning: Unknown option "$1" + echo rerun with --help for more info. + exit 1 + fi + + shift 1 + +done + +if [ $# -eq 1 ]; then + echo "Prepending $1 to your PATH for the duration of this script." + PATH=$1:$PATH +fi + +echo "Using `command -v SSRB`" +echo "Using `command -v OSMAPOSL`" + +# first need to set this to the C locale, as this is what the STIR utilities use +# otherwise, awk might interpret floating point numbers incorrectly +LC_ALL=C +export LC_ALL + +echo "=== Creating single view data and test TOF rebinning ===" +logfile=SSRB_sim1.1og +max_rd=9 span=3 view_mash=288 tof_mash=1 ./simulate_PET_data_for_tests.sh --suffix _TOF_vm288_rd9_span3 --TOF > "$logfile" 2>&1 +if [ $? -ne 0 ]; then + echo "Error running simulation, Check $logfile" + exit 1 +fi +logfile=SSRB_sim2.1og +max_rd=9 span=3 view_mash=288 tof_mash=11 ./simulate_PET_data_for_tests.sh --keep_images --suffix _TOF_vm288_rd9_span3_tm11 --TOF > "$logfile" 2>&1 +if [ $? -ne 0 ]; then + echo "Error running simulation, Check $logfile" + exit 1 +fi +logfile=SSRB_run_SSRB_sim1.log +SSRB --template my_line_integrals_TOF_vm288_rd9_span3_tm11.hs my_line_integrals_TOF_vm288_rd9_span3_TOFSSRB11.hs my_line_integrals_TOF_vm288_rd9_span3.hs 0 > "$logfile" 2>&1 +if [ $? -ne 0 ]; then + echo "Error running SSRB, Check $logfile" + exit 1 +fi +compare_projdata my_line_integrals_TOF_vm288_rd9_span3_tm11.hs my_line_integrals_TOF_vm288_rd9_span3_TOFSSRB11.hs +if [ $? -ne 0 ]; then + echo "Error comparing TOF mash simulation with TOF-SSRB" + exit 1 +fi + + +echo "" +echo "=== Creating 3D data and test axial+TOF rebinning via reconstruction ===" +echo "We will simulate some data at span=1 and TOF-mashing=1, rebin it, then reconstruct." +echo "The image will be checked by using an ROI calculation." + +echo "Running simulation (will take a while)" +max_rd=5 +export max_rd +view_mash=8 +export view_mash +suffix="_TOF_vm${view_mash}_rd${max_rd}_span1" +logfile=SSRB_sim3.1og +span=1 tof_mash=1 background_value=.1 ./simulate_PET_data_for_tests.sh --keep_images --suffix "$suffix" --TOF > "$logfile" 2>&1 +if [ $? -ne 0 ]; then + echo "Error running simulation, Check $logfile" + exit 1 +fi +echo "Rebinning data" +tof_mash=11 +span=3 +suffixSSRB="_TOF_vm${view_mash}_rd${max_rd}_span${span}_TOFSSRB${tof_mash}" +logfile="SSRB_run_SSRB_line_integrals$suffixSSRB.1og" +SSRB "my_line_integrals${suffixSSRB}.hs" "my_line_integrals${suffix}.hs" "$span" 1 1 "$max_rd" "$tof_mash" > "$logfile" 2>&1 +if [ $? -ne 0 ]; then + echo "Error running SSRB, Check $logfile" + exit 1 +fi +logfile="SSRB_run_SSRB_acfs$suffixSSRB.1og" +SSRB "my_acfs${suffixSSRB}.hs" "my_acfs${suffix}.hs" "$span" 1 1 "$max_rd" "$tof_mash" > "$logfile" 2>&1 +if [ $? -ne 0 ]; then + echo "Error running SSRB, Check $logfile" + exit 1 +fi +logfile="SSRB_run_SSRB_prompts$suffixSSRB.1og" +SSRB "my_prompts${suffixSSRB}.hs" "my_prompts${suffix}.hs" "$span" 1 1 "$max_rd" "$tof_mash" > "$logfile" 2>&1 +if [ $? -ne 0 ]; then + echo "Error running SSRB, Check $logfile" + exit 1 +fi +logfile="SSRB_run_SSRB_additive_sinograms$suffixSSRB.1og" +SSRB "my_additive_sinogram${suffixSSRB}.hs" "my_additive_sinogram${suffix}.hs" "$span" 1 1 "$max_rd" "$tof_mash" > "$logfile" 2>&1 +if [ $? -ne 0 ]; then + echo "Error running SSRB, Check $logfile" + exit 1 +fi + + +# Below code is a copy of relevant lines in run_test_simulate_and_recon.sh. +# It would be far better to make those into functions/scripts... + +recon=OSMAPOSL +isFBP=0 +parfile=${recon}_test_sim_PM.par +suffix=${suffixSSRB} +export suffix + +num_subsets=2 +export num_subsets + +input_image=my_uniform_cylinder.hv +input_voxel_size_x=`stir_print_voxel_sizes.sh ${input_image}|awk '{print $3}'` +ROI=ROI_uniform_cylinder.par +list_ROI_values ${input_image}.roistats ${input_image} ${ROI} 0 > /dev/null 2>&1 +input_ROI_mean=`awk 'NR>2 {print $2}' ${input_image}.roistats` + +# Indentation is here because straight copy of the file. + + # run actual reconstruction + echo "Running ${recon} ${parfile}" + logfile="my_${parfile}${suffix}.log" + ${MPIRUN} ${recon} ${parfile} > "$logfile" 2>&1 + if [ $? -ne 0 ]; then + echo "Error running reconstruction. CHECK RECONSTRUCTION LOG \"$logfile\"" + exit 1 + fi + + # find filename of (last) image from ${parfile} + output_filename=`awk -F':=' '/output[ _]*filename[ _]*prefix/ { value=$2;gsub(/[ \t]/, "", value); printf("%s", value) }' "$parfile"` + # substitute env variables (e.g. to fill in suffix) + output_filename=`eval echo "${output_filename}"` + if [ ${isFBP} -eq 0 ]; then + # iterative algorithm, so we need to append the num_subiterations + num_subiterations=`awk -F':=' '/number[ _]*of[ _]*subiterations/ { value=$2;gsub(/[ \t]/, "", value); printf("%s", value) }' ${parfile}` + output_filename=${output_filename}_${num_subiterations} + fi + output_image=${output_filename}.hv + + # compute ROI value + list_ROI_values ${output_image}.roistats ${output_image} ${ROI} 0 > ${output_image}.roistats.log 2>&1 + if [ $? -ne 0 ]; then + echo "Error running list_ROI_values. CHECK LOG ${output_image}.roistats.log" + exit 1 + fi + + # compare ROI value + output_voxel_size_x=`stir_print_voxel_sizes.sh ${output_image}|awk '{print $3}'` + output_ROI_mean=`awk "NR>2 {print \\$2*${input_voxel_size_x}/${output_voxel_size_x}}" ${output_image}.roistats` + echo "Input ROI mean: $input_ROI_mean" + echo "Output ROI mean: $output_ROI_mean" + error_bigger_than_1percent=`echo $input_ROI_mean $output_ROI_mean| awk '{ print(($2/$1 - 1)*($2/$1 - 1)>0.0001) }'` + if [ ${error_bigger_than_1percent} -eq 1 ]; then + echo "DIFFERENCE IN ROI VALUES IS TOO LARGE. CHECK RECONSTRUCTION LOG "$logfile"" + else + echo "This seems fine." + fi + From 24a910259fdb331834e8aaccacf52c435330f8da Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 4 Jul 2024 13:26:20 +0100 Subject: [PATCH 16/17] convert some old-style warning() to new-style When using a single int as argument, the old style version is actually not used anymore due to the overload with the verbosity presumably. So use boost::format --- src/buildblock/ProjDataInfoCylindrical.cxx | 31 +++++++++++----------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/src/buildblock/ProjDataInfoCylindrical.cxx b/src/buildblock/ProjDataInfoCylindrical.cxx index 8cc5942cbc..e9e4e29566 100644 --- a/src/buildblock/ProjDataInfoCylindrical.cxx +++ b/src/buildblock/ProjDataInfoCylindrical.cxx @@ -107,11 +107,10 @@ ProjDataInfoCylindrical::ProjDataInfoCylindrical(const shared_ptr& scan for (int segment_num = get_min_segment_num(); segment_num <= get_max_segment_num(); ++segment_num) if (min_ring_diff[segment_num] > max_ring_diff[segment_num]) { - warning("ProjDataInfoCylindrical: min_ring_difference %d is larger than max_ring_difference %d for segment %d. " - "Swapping them around", - min_ring_diff[segment_num], - max_ring_diff[segment_num], - segment_num); + warning(boost::format( + "ProjDataInfoCylindrical: min_ring_difference %d is larger than max_ring_difference %d for segment %d. " + "Swapping them around") + % min_ring_diff[segment_num] % max_ring_diff[segment_num] % segment_num); std::swap(min_ring_diff[segment_num], max_ring_diff[segment_num]); } } @@ -199,11 +198,11 @@ ProjDataInfoCylindrical::initialise_ring_diff_arrays() const // check that it was integer if (fabs(ax_pos_num_offset[segment_num] - ((num_rings - 1) - 2 * m_offset[segment_num] / ring_spacing)) > 1E-4) { - error("ProjDataInfoCylindrical: in segment %d, the axial positions\n" - "do not correspond to the usual locations between physical rings.\n" - "This is suspicious and can make things go wrong in STIR, so I abort.\n" - "Check the number of axial positions in this segment.", - segment_num); + error(boost::format("ProjDataInfoCylindrical: in segment %d, the axial positions\n" + "do not correspond to the usual locations between physical rings.\n" + "This is suspicious and can make things go wrong in STIR, so I abort.\n" + "Check the number of axial positions in this segment.") + % segment_num); } } @@ -215,10 +214,12 @@ ProjDataInfoCylindrical::initialise_ring_diff_arrays() const // ring1+ring2 = 2*ring2 + ring_diff assert(get_min_ring_difference(segment_num) == get_max_ring_difference(segment_num)); if ((get_max_ring_difference(segment_num) - ax_pos_num_offset[segment_num]) % 2 != 0) - warning("ProjDataInfoCylindrical: the number of axial positions in " - "segment %d is such that current conventions will place " - "the LORs shifted with respect to the physical rings.", - segment_num); + warning(boost::format("ProjDataInfoCylindrical: the number of axial positions %d in " + "segment %d (ring_diff %d) is such that current conventions will place " + "the LORs shifted with respect to the physical rings %d.") + % get_num_axial_poss(segment_num) % segment_num + % get_min_ring_difference(segment_num) // equal to max here, as per the if() + % get_scanner_ptr()->get_num_rings()); } } } @@ -255,7 +256,7 @@ ProjDataInfoCylindrical::initialise_ring_diff_arrays() const } if (segment_num > get_max_segment_num()) { - warning("ProjDataInfoCylindrical: ring difference %d does not belong to a segment", ring_diff); + warning(boost::format("ProjDataInfoCylindrical: ring difference %d does not belong to a segment") % ring_diff); } } } From 131f3ad598d2c05f7dcbf7f28285aa364ea4d955 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 4 Jul 2024 13:27:41 +0100 Subject: [PATCH 17/17] change max_rd for STE test-case span=2, max_rd=2 is unexpected, and leads to a bug in FBP3DRP etc. Changing to the expected case of max_rd=3 for now. --- recon_test_pack/simulate_PET_data_for_tests.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/recon_test_pack/simulate_PET_data_for_tests.sh b/recon_test_pack/simulate_PET_data_for_tests.sh index 94e884a0bf..6f31b26189 100755 --- a/recon_test_pack/simulate_PET_data_for_tests.sh +++ b/recon_test_pack/simulate_PET_data_for_tests.sh @@ -91,7 +91,7 @@ fi if [ "$TOF" -eq 0 ]; then : ${view_mash:=1} : ${span:=2} - : ${max_rd:=2} + : ${max_rd:=3} echo "=== create template sinogram (DSTE with view_mash=${view_mash}, max_ring_diff=${max_rd})" template_sino=my_DSTE_3D_vm${view_mash}_span${span}_rd${max_rd}_template.hs cat > my_input.txt <