From 59b812db44cb26b93f07f6b620d7d9627d3e9dd9 Mon Sep 17 00:00:00 2001 From: Somayeh Date: Wed, 14 Feb 2024 10:48:49 +0100 Subject: [PATCH 01/19] Initial commit: added Safir-I/II to Scanner.cxx/.h --- src/buildblock/Scanner.cxx | 58 ++++++++++++++++++++++++++++++++++++++ src/include/stir/Scanner.h | 2 ++ 2 files changed, 60 insertions(+) diff --git a/src/buildblock/Scanner.cxx b/src/buildblock/Scanner.cxx index a66b5c3d4b..1d78e8a3c1 100644 --- a/src/buildblock/Scanner.cxx +++ b/src/buildblock/Scanner.cxx @@ -1110,6 +1110,64 @@ Scanner::Scanner(Type scanner_type) ); break; + case SafirI: + set_params(SafirI, string_list("SafirI"), + 24, //num_rings_v + 150, //max_num_non_arccorrected_bins_v, + 150, //default_num_arccorrected_bins_v, + 180, //num_detectors_per_ring_v + 64.05, // inner_ring_radius_v + 5, //average_depth_of_interaction_v + 2.2, //ring_spacing_v + 1.1, //bin_size_v + 0, //intrinsic_tilt_v + 3, //num_axial_blocks_per_bucket_v + 1, //num_transaxial_blocks_per_bucket_v + 8, //num_axial_crystals_per_block_v + 15, //num_transaxial_crystals_per_block_v + 1, //num_axial_crystals_per_singles_unit_v + 1, //num_transaxial_crystals_per_singles_unit_v + 1, //num_detector_layers_v + 0.11, //energy_resolution_v + 511, //reference_energy_v + "", //scanner_geometry_v + 2.2, //axial_crystal_spacing_v + 2.2, //transaxial_crystal_spacing_v + 18.1, //axial_block_spacing_v + 33.6, //transaxial_block_spacing_v + ""//crystal_map_file_name_v + ); + break; + + case SafirII: + set_params(SafirII, string_list("SafirII"), + 64, //num_rings_v + 150, //max_num_non_arccorrected_bins_v, + 150, //default_num_arccorrected_bins_v, + 180, //num_detectors_per_ring_v + 64.05, // inner_ring_radius_v + 5, //average_depth_of_interaction_v + 2.2, //ring_spacing_v + 1.1, //bin_size_v + 0, //intrinsic_tilt_v + 8, //num_axial_blocks_per_bucket_v + 1, //num_transaxial_blocks_per_bucket_v + 8, //num_axial_crystals_per_block_v + 15, //num_transaxial_crystals_per_block_v + 1, //num_axial_crystals_per_singles_unit_v + 1, //num_transaxial_crystals_per_singles_unit_v + 1, //num_detector_layers_v + 0.11, //energy_resolution_v + 511, //reference_energy_v + "", //scanner_geometry_v + 2.2, //axial_crystal_spacing_v + 2.2, //transaxial_crystal_spacing_v + 18.1, //axial_block_spacing_v + 33.6, //transaxial_block_spacing_v + ""//crystal_map_file_name_v + ); + break; + case UPENN_5rings: set_params(UPENN_5rings, string_list("UPENN_5rings"), diff --git a/src/include/stir/Scanner.h b/src/include/stir/Scanner.h index 79c8cbde82..a1b1b40ba6 100644 --- a/src/include/stir/Scanner.h +++ b/src/include/stir/Scanner.h @@ -174,6 +174,8 @@ class Scanner Allegro, GeminiTF, SAFIRDualRingPrototype, + SafirI, + SafirII, UPENN_5rings, UPENN_5rings_no_gaps, UPENN_6rings, From fc76b4060ff198a372ceeb5801427f4be08874f4 Mon Sep 17 00:00:00 2001 From: Somayeh Date: Wed, 14 Feb 2024 10:49:52 +0100 Subject: [PATCH 02/19] added a new branch build folder to gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index d3689c4f86..b183a3bebd 100644 --- a/.gitignore +++ b/.gitignore @@ -70,6 +70,7 @@ install_manifest.txt build/ install/ .vscode/ +build_safir_updates/ # MATLAB From 4407e3b396e0d7e5a0fe345d89e06ef20048a5f4 Mon Sep 17 00:00:00 2001 From: ssaghamanesh Date: Thu, 5 Sep 2024 11:20:26 +0200 Subject: [PATCH 03/19] Added Safir I/II scanners --- src/buildblock/Scanner.cxx | 84 +++++++++++++++++++++++++++++++++----- src/include/stir/Scanner.h | 2 + 2 files changed, 76 insertions(+), 10 deletions(-) diff --git a/src/buildblock/Scanner.cxx b/src/buildblock/Scanner.cxx index 2f48f14387..9ddcab555c 100644 --- a/src/buildblock/Scanner.cxx +++ b/src/buildblock/Scanner.cxx @@ -1085,7 +1085,7 @@ Scanner::Scanner(Type scanner_type) 150, // max_num_non_arccorrected_bins_v, 150, // default_num_arccorrected_bins_v, 180, // num_detectors_per_ring_v - 64.05, // inner_ring_radius_v + 64.05, // inner_ring_radius_v 5, // average_depth_of_interaction_v 2.2, // ring_spacing_v 1.1, // bin_size_v @@ -1099,18 +1099,82 @@ Scanner::Scanner(Type scanner_type) 1, // num_detector_layers_v -1, // energy_resolution_v -1, // reference_energy_v - (short int)1, - 0.F, - 0.F, // non-TOF - "BlocksOnCylindrical", // scanner_geometry_v - 2.2, // axial_crystal_spacing_v - 2.2, // transaxial_crystal_spacing_v - 18.1, // axial_block_spacing_v - 33.6, // transaxial_block_spacing_v - "" // crystal_map_file_name_v + (short int)1, // max_num_of_timing_poss_v, + 0.F, // size_timing_pos_v, + 0.F, // timing_resolution_v, + "", // scanner_geometry_v + 2.2, // axial_crystal_spacing_v + 2.2, // transaxial_crystal_spacing_v + 18.1, // axial_block_spacing_v + 33.6, // transaxial_block_spacing_v + "" // crystal_map_file_name_v ); break; + case SafirI: + set_params(SafirI, string_list("SafirI"), + 24, // num_rings_v + 150, // max_num_non_arccorrected_bins_v, + 150, // default_num_arccorrected_bins_v, + 180, // num_detectors_per_ring_v + 64.05, // inner_ring_radius_v + 5, // average_depth_of_interaction_v + 2.2, // ring_spacing_v + 1.1, // bin_size_v + 0, // intrinsic_tilt_v + 3, // num_axial_blocks_per_bucket_v + 1, // num_transaxial_blocks_per_bucket_v + 8, // num_axial_crystals_per_block_v + 15, // num_transaxial_crystals_per_block_v + 1, // num_axial_crystals_per_singles_unit_v + 1, // num_transaxial_crystals_per_singles_unit_v + 1, // num_detector_layers_v + 0.12, // energy_resolution_v + 511, // reference_energy_v + (short int)1, //max_num_of_timing_poss_v, + 0.F, // size_timing_pos_v, + 0.F, // timing_resolution_v, + "", // scanner_geometry_v + 2.2, // axial_crystal_spacing_v + 2.2, // transaxial_crystal_spacing_v + 18.1, // axial_block_spacing_v + 33.6, // transaxial_block_spacing_v + "" // crystal_map_file_name_v + ); + break; + + case SafirII: + set_params(SafirII, string_list("SafirII"), + 64, // num_rings_v + 150, // max_num_non_arccorrected_bins_v, + 150, // default_num_arccorrected_bins_v, + 180, // num_detectors_per_ring_v + 64.05, // inner_ring_radius_v + 5, // average_depth_of_interaction_v + 2.2, // ring_spacing_v + 1.1, // bin_size_v + 0, // intrinsic_tilt_v + 8, // num_axial_blocks_per_bucket_v + 1, // num_transaxial_blocks_per_bucket_v + 8, // num_axial_crystals_per_block_v + 15, // num_transaxial_crystals_per_block_v + 1, // num_axial_crystals_per_singles_unit_v + 1, // num_transaxial_crystals_per_singles_unit_v + 1, // num_detector_layers_v + 0.12, // energy_resolution_v + 511, // reference_energy_v + (short int)1, // max_num_of_timing_poss_v, + 0.F, // size_timing_pos_v, + 0.F, // timing_resolution_v, + "", // scanner_geometry_v + 2.2, // axial_crystal_spacing_v + 2.2, // transaxial_crystal_spacing_v + 18.1, // axial_block_spacing_v + 33.6, // transaxial_block_spacing_v + "" // crystal_map_file_name_v + ); + break; + case UPENN_5rings: set_params(UPENN_5rings, string_list("UPENN_5rings"), diff --git a/src/include/stir/Scanner.h b/src/include/stir/Scanner.h index e5b7332bd7..97ee665a3d 100644 --- a/src/include/stir/Scanner.h +++ b/src/include/stir/Scanner.h @@ -174,6 +174,8 @@ class Scanner Allegro, GeminiTF, SAFIRDualRingPrototype, + SafirI, + SafirII, UPENN_5rings, UPENN_5rings_no_gaps, UPENN_6rings, From 69486f5ac29c41a9ece446f6b0efda7ee52967ed Mon Sep 17 00:00:00 2001 From: ssaghamanesh Date: Thu, 5 Sep 2024 11:21:07 +0200 Subject: [PATCH 04/19] Adding two utilities for conversion and trim sinograms in generic geometry --- src/utilities/convert_projdata_types.cxx | 95 +++++++++ src/utilities/trim_projdata.cxx | 246 +++++++++++++++++++++++ 2 files changed, 341 insertions(+) create mode 100644 src/utilities/convert_projdata_types.cxx create mode 100644 src/utilities/trim_projdata.cxx diff --git a/src/utilities/convert_projdata_types.cxx b/src/utilities/convert_projdata_types.cxx new file mode 100644 index 0000000000..dbbac8f157 --- /dev/null +++ b/src/utilities/convert_projdata_types.cxx @@ -0,0 +1,95 @@ +// +// +/*! + \file + \ingroup utilities + + \brief This program takes a projection data from one type and converts it to another type. + \author Parisa Khateri + +*/ +/* +*/ +#include "stir/ProjData.h" +#include "stir/IO/interfile.h" +#include "stir/utilities.h" +#include "stir/Bin.h" + +#include +#include +#include "stir/ProjDataFromStream.h" +#include "stir/Viewgram.h" +#include "stir/IO/read_from_file.h" +#include "stir/SegmentByView.h" +#include "stir/ProjDataInterfile.h" +#include "stir/ProjDataInfo.h" +#include "stir/LORCoordinates.h" + +#include "stir/GeometryBlocksOnCylindrical.h" +#include "stir/DetectionPosition.h" +#include "stir/CartesianCoordinate3D.h" +#include "stir/listmode/DetectorCoordinateMapFromFile.h" +#include +#include "stir/CPUTimer.h" +#include "stir/shared_ptr.h" + + +#ifndef STIR_NO_NAMESPACES +using std::cerr; +#endif + +USING_NAMESPACE_STIR + + + +int main(int argc, char *argv[]) +{ + CPUTimer timer0; + + if(argc<4) + { + cerr<<"Usage: " << argv[0] << " output_filename input_filename template_blk_projdata\n"; + exit(EXIT_FAILURE); + } + std::string output_filename=argv[1]; + shared_ptr in_pd_ptr = ProjData::read_from_file(argv[2]); + shared_ptr template_pd_ptr = ProjData::read_from_file(argv[3]); + + shared_ptr in_pdi_ptr(in_pd_ptr->get_proj_data_info_sptr()->clone()); + shared_ptr out_pdi_ptr(template_pd_ptr->get_proj_data_info_sptr()->clone()); + ProjDataInterfile out_proj_data(template_pd_ptr->get_exam_info_sptr(), out_pdi_ptr, output_filename+".hs"); + write_basic_interfile_PDFS_header(output_filename+".hs", out_proj_data); + + timer0.start(); + + assert(in_pdi_ptr->get_min_segment_num()==-1*in_pdi_ptr->get_max_segment_num()); + for (int seg=in_pdi_ptr->get_min_segment_num(); seg<=in_pdi_ptr->get_max_segment_num();++seg) + { + std::cout<<"seg_num = "< viewgram_blk = out_proj_data.get_empty_viewgram(out_proj_data.get_min_view_num(),seg); + Viewgram viewgram_cyl = in_pd_ptr->get_empty_viewgram(in_pd_ptr->get_min_view_num(),seg); + + for(int view=in_pdi_ptr->get_min_view_num(); view<=in_pdi_ptr->get_max_view_num();++view) + { + viewgram_blk = out_proj_data.get_empty_viewgram(view,seg); + viewgram_cyl = in_pd_ptr->get_viewgram(view,seg); + + for(int ax=in_pdi_ptr->get_min_axial_pos_num(seg); ax<=in_pdi_ptr->get_max_axial_pos_num(seg);++ax) + { + for(int tang=in_pdi_ptr->get_min_tangential_pos_num(); tang<=in_pdi_ptr->get_max_tangential_pos_num(); ++tang) + { + viewgram_blk[ax][tang] = viewgram_cyl[ax][tang]; + } + } + if (!(out_proj_data.set_viewgram(viewgram_blk)== Succeeded::yes)) + warning("Error set_segment for projdata_symm %d\n", seg); + } + } + + timer0.stop(); + std::cerr << "\nConverting from cylindrical projdata to block took " << timer0.value() << "s CPU time.\n\n"; + + return EXIT_SUCCESS; +} diff --git a/src/utilities/trim_projdata.cxx b/src/utilities/trim_projdata.cxx new file mode 100644 index 0000000000..208f899046 --- /dev/null +++ b/src/utilities/trim_projdata.cxx @@ -0,0 +1,246 @@ +// +// +/* + +*/ +/*! + + \file + \ingroup utilities + \brief Main program for trim_projdata + + \author Parisa Khateri + + + \par Usage: + \code + trim_projdata [-t num_tang_poss_to_trim] \ + output_filename input_projdata_name + \endcode + \param num_tang_poss_to_trim has to be smaller than the available number + of tangential positions. + + \par Example: + \code + trim_projdata -t 10 out in.hs + \endcode +*/ +#include "stir/ProjData.h" +#include "stir/shared_ptr.h" +#include +#include +#include "stir/ProjDataFromStream.h" +#include "stir/ProjDataInterfile.h" +#include "stir/ProjDataInfoCylindrical.h" +#include "stir/ProjDataInfoBlocksOnCylindrical.h" +#include "stir/ProjDataInfoGeneric.h" +#include "stir/Sinogram.h" +#include "stir/Bin.h" +#include "stir/round.h" +#include +#include + + +#ifndef STIR_NO_NAMESPACES +using std::string; +using std::cerr; +#endif + +USING_NAMESPACE_STIR + +int main(int argc, char **argv) +{ + int num_tang_poss_to_trim = 0; + if (argc>1 && strcmp(argv[1], "-t")==0) + { + num_tang_poss_to_trim = atoi(argv[2]); + argc -= 2; argv += 2; + } + if (argc > 5 || argc < 3 ) + { + cerr << "Usage:\n" + << argv[0] << " [-t num_tang_poss_to_trim] \\\n" + << "\toutput_filename input_projdata_name \\\n" + << "num_tang_poss_to_trim has to be smaller than the available number\n"; + exit(EXIT_FAILURE); + } + const string output_filename = argv[1]; + shared_ptr in_projdata_ptr = ProjData::read_from_file(argv[2]); + + + if (in_projdata_ptr->get_num_tangential_poss() <= + num_tang_poss_to_trim) + error("trim_projdata: too large number of tangential positions to trim (%d)\n", + num_tang_poss_to_trim); + + if (in_projdata_ptr->get_proj_data_info_sptr()->get_scanner_ptr()->get_scanner_geometry() == + "Cylindrical") + { + const ProjDataInfoCylindrical * const in_projdata_info_cyl_ptr = + dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + if (in_projdata_info_cyl_ptr== NULL) + { + error("error converting to cylindrical projection data\n"); + } + ProjDataInfoCylindrical * out_projdata_info_cyl_ptr = + dynamic_cast + (in_projdata_info_cyl_ptr->clone()); + + out_projdata_info_cyl_ptr-> + set_num_tangential_poss(in_projdata_info_cyl_ptr->get_num_tangential_poss() - + num_tang_poss_to_trim); + + shared_ptr out_projdata_info_ptr(out_projdata_info_cyl_ptr); + ProjDataInterfile out_projdata(in_projdata_ptr->get_exam_info_sptr(), + out_projdata_info_ptr, output_filename, std::ios::out); + + for (int seg = out_projdata.get_min_segment_num(); + seg <= out_projdata.get_max_segment_num(); + ++seg) + { + // keep sinograms out of the loop to avoid reallocations + // initialise to something because there's no default constructor + Sinogram out_sino = + out_projdata.get_empty_sinogram(out_projdata.get_min_axial_pos_num(seg),seg); + Sinogram in_sino = + in_projdata_ptr->get_empty_sinogram(in_projdata_ptr->get_min_axial_pos_num(seg),seg); + + for (int ax = out_projdata.get_min_axial_pos_num(seg); + ax <= out_projdata.get_max_axial_pos_num(seg); + ++ax ) + { + out_sino= out_projdata.get_empty_sinogram(ax, seg); + in_sino = in_projdata_ptr->get_sinogram(ax, seg); + + { + for (int view=out_projdata.get_min_view_num(); + view <= out_projdata.get_max_view_num(); + ++view) + for (int tang=out_projdata.get_min_tangential_pos_num(); + tang <= out_projdata.get_max_tangential_pos_num(); + ++tang) + out_sino[view][tang] = in_sino[view][tang]; + } + out_projdata.set_sinogram(out_sino); + } + + } + + } + else if (in_projdata_ptr->get_proj_data_info_sptr()->get_scanner_ptr()->get_scanner_geometry() == + "BlocksOnCylindrical") + { + const ProjDataInfoBlocksOnCylindrical * const in_projdata_info_blk_ptr = + dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + if (in_projdata_info_blk_ptr== NULL) + { + error("error converting to BlocksOnCylindrical projection data\n"); + } + ProjDataInfoBlocksOnCylindrical * out_projdata_info_blk_ptr = + dynamic_cast + (in_projdata_info_blk_ptr->clone()); + + out_projdata_info_blk_ptr-> + set_num_tangential_poss(in_projdata_info_blk_ptr->get_num_tangential_poss() - + num_tang_poss_to_trim); + + shared_ptr out_projdata_info_ptr(out_projdata_info_blk_ptr); + ProjDataInterfile out_projdata(in_projdata_ptr->get_exam_info_sptr(), + out_projdata_info_ptr, output_filename, std::ios::out); + + + for (int seg = out_projdata.get_min_segment_num(); + seg <= out_projdata.get_max_segment_num(); + ++seg) + { + // keep sinograms out of the loop to avoid reallocations + // initialise to something because there's no default constructor + Sinogram out_sino = + out_projdata.get_empty_sinogram(out_projdata.get_min_axial_pos_num(seg),seg); + Sinogram in_sino = + in_projdata_ptr->get_empty_sinogram(in_projdata_ptr->get_min_axial_pos_num(seg),seg); + + for (int ax = out_projdata.get_min_axial_pos_num(seg); + ax <= out_projdata.get_max_axial_pos_num(seg); + ++ax ) + { + out_sino= out_projdata.get_empty_sinogram(ax, seg); + in_sino = in_projdata_ptr->get_sinogram(ax, seg); + + { + for (int view=out_projdata.get_min_view_num(); + view <= out_projdata.get_max_view_num(); + ++view) + for (int tang=out_projdata.get_min_tangential_pos_num(); + tang <= out_projdata.get_max_tangential_pos_num(); + ++tang) + out_sino[view][tang] = in_sino[view][tang]; + } + out_projdata.set_sinogram(out_sino); + } + + } + + + } + else if (in_projdata_ptr->get_proj_data_info_sptr()->get_scanner_ptr()->get_scanner_geometry() == + "Generic") + { + const ProjDataInfoGeneric * const in_projdata_info_gen_ptr = + dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + if (in_projdata_info_gen_ptr== NULL) + { + error("error converting to Generic projection data\n"); + } + ProjDataInfoGeneric * out_projdata_info_gen_ptr = + dynamic_cast + (in_projdata_info_gen_ptr->clone()); + + out_projdata_info_gen_ptr-> + set_num_tangential_poss(in_projdata_info_gen_ptr->get_num_tangential_poss() - + num_tang_poss_to_trim); + + shared_ptr out_projdata_info_ptr(out_projdata_info_gen_ptr); + ProjDataInterfile out_projdata(in_projdata_ptr->get_exam_info_sptr(), + out_projdata_info_ptr, output_filename, std::ios::out); + + for (int seg = out_projdata.get_min_segment_num(); + seg <= out_projdata.get_max_segment_num(); + ++seg) + { + // keep sinograms out of the loop to avoid reallocations + // initialise to something because there's no default constructor + Sinogram out_sino = + out_projdata.get_empty_sinogram(out_projdata.get_min_axial_pos_num(seg),seg); + Sinogram in_sino = + in_projdata_ptr->get_empty_sinogram(in_projdata_ptr->get_min_axial_pos_num(seg),seg); + + for (int ax = out_projdata.get_min_axial_pos_num(seg); + ax <= out_projdata.get_max_axial_pos_num(seg); + ++ax ) + { + out_sino= out_projdata.get_empty_sinogram(ax, seg); + in_sino = in_projdata_ptr->get_sinogram(ax, seg); + + { + for (int view=out_projdata.get_min_view_num(); + view <= out_projdata.get_max_view_num(); + ++view) + for (int tang=out_projdata.get_min_tangential_pos_num(); + tang <= out_projdata.get_max_tangential_pos_num(); + ++tang) + out_sino[view][tang] = in_sino[view][tang]; + } + out_projdata.set_sinogram(out_sino); + } + + } + + } + else + { + error("error the scanner geometry of projection data is not known\n"); + } + + return EXIT_SUCCESS; +} From bf952407699c8a9373e6183064fedd304a4e5fce Mon Sep 17 00:00:00 2001 From: ssaghamanesh Date: Thu, 5 Sep 2024 15:12:11 +0200 Subject: [PATCH 05/19] Excluded .gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 075f775d2a..400894565f 100644 --- a/.gitignore +++ b/.gitignore @@ -71,6 +71,8 @@ install_manifest.txt build/ install/ .vscode/ +build_safir_updates/ +.gitignore # MATLAB From 65802345de34a0beba3685068111ffa44e67825b Mon Sep 17 00:00:00 2001 From: ssaghamanesh Date: Thu, 5 Sep 2024 15:27:32 +0200 Subject: [PATCH 06/19] Initial commit + excluded build in .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 075f775d2a..51c7038fde 100644 --- a/.gitignore +++ b/.gitignore @@ -71,6 +71,7 @@ install_manifest.txt build/ install/ .vscode/ +build_safir_updates/ # MATLAB From 566100bfe8c5ba1e4e7e68c2aec2b4ef09d1153a Mon Sep 17 00:00:00 2001 From: ssaghamanesh Date: Thu, 5 Sep 2024 15:28:18 +0200 Subject: [PATCH 07/19] Modified min/max_tang_pos_num to be consistent with ProjDataInfo.cxx --- src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx | 4 ++-- src/buildblock/ProjDataInfoGenericNoArcCorr.cxx | 4 ++-- src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx | 4 +++- 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx index 93b2e49169..a11543cd7a 100644 --- a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx @@ -178,8 +178,8 @@ ProjDataInfoCylindricalNoArcCorr::initialise_uncompressed_view_tangpos_to_det1de assert(fabs(get_phi(Bin(0, 0, 0, 0)) - v_offset) < 1.E-4); assert(fabs(get_phi(Bin(0, get_num_views(), 0, 0)) - v_offset - _PI) < 1.E-4); #endif - const int min_tang_pos_num = -(num_detectors / 2) + 1; - const int max_tang_pos_num = -(num_detectors / 2) + num_detectors; + const int min_tang_pos_num = -(num_detectors / 2); + const int max_tang_pos_num = -(num_detectors / 2) + num_detectors - 1; if (this->get_min_tangential_pos_num() < min_tang_pos_num || this->get_max_tangential_pos_num() > max_tang_pos_num) { diff --git a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx index 6b4ea7de48..930e3be73a 100644 --- a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx @@ -149,8 +149,8 @@ ProjDataInfoGenericNoArcCorr::initialise_uncompressed_view_tangpos_to_det1det2() const int num_detectors = get_scanner_ptr()->get_num_detectors_per_ring(); assert(num_detectors % 2 == 0); - const int min_tang_pos_num = -(num_detectors / 2) + 1; - const int max_tang_pos_num = -(num_detectors / 2) + num_detectors; + const int min_tang_pos_num = -(num_detectors / 2); + const int max_tang_pos_num = -(num_detectors / 2) + num_detectors - 1; if (this->get_min_tangential_pos_num() < min_tang_pos_num || this->get_max_tangential_pos_num() > max_tang_pos_num) { diff --git a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx index 3d596800bb..5613e3fe57 100644 --- a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx +++ b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx @@ -273,9 +273,11 @@ ProjMatrixByBinUsingRayTracing::set_up( voxel_size = image_info_ptr->get_voxel_size(); origin = image_info_ptr->get_origin(); - if (std::abs(origin.x()) > .05F || std::abs(origin.y()) > .05F) + if (std::abs(origin.x()) > .05F || std::abs(origin.y()) > .05F){ error("ProjMatrixByBinUsingRayTracing sadly doesn't support shifted x/y origin yet"); image_info_sptr->get_regular_range(min_index, max_index); + std::cout << "OriginX : " << origin.x() << " ; OriginY : " << origin.y() << "\n"; + } symmetries_sptr.reset(new DataSymmetriesForBins_PET_CartesianGrid(proj_data_info_sptr, density_info_sptr_v, From 8ae5dd5f116491f9b055149166bdf0c718adacd4 Mon Sep 17 00:00:00 2001 From: Somayeh Saghamanesh <96595217+SomayehSaghamanesh@users.noreply.github.com> Date: Fri, 6 Sep 2024 17:33:22 +0200 Subject: [PATCH 08/19] Updated src/utilities/CMakeLists.txt --- src/utilities/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt index 42e09b8318..1d99a99dfa 100644 --- a/src/utilities/CMakeLists.txt +++ b/src/utilities/CMakeLists.txt @@ -75,6 +75,8 @@ set(${dir_EXE_SOURCES} find_sum_projection_of_viewgram_and_sinogram.cxx separate_true_from_random_scatter_for_necr.cxx stir_timings.cxx + convert_projdata_types.cxx + trim_projdata.cxx ) if (HAVE_ITK) From 5f4d9355d3d8d8181b52dcc5d6d62338505ca9ce Mon Sep 17 00:00:00 2001 From: ssaghamanesh Date: Mon, 9 Sep 2024 12:13:04 +0200 Subject: [PATCH 09/19] Fixed dynamic_cast problem in two utilities --- src/utilities/convert_projdata_types.cxx | 2 +- src/utilities/trim_projdata.cxx | 19 +++++++++++++------ 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/src/utilities/convert_projdata_types.cxx b/src/utilities/convert_projdata_types.cxx index dbbac8f157..d254da003c 100644 --- a/src/utilities/convert_projdata_types.cxx +++ b/src/utilities/convert_projdata_types.cxx @@ -28,7 +28,7 @@ #include "stir/GeometryBlocksOnCylindrical.h" #include "stir/DetectionPosition.h" #include "stir/CartesianCoordinate3D.h" -#include "stir/listmode/DetectorCoordinateMapFromFile.h" +#include "stir/DetectorCoordinateMap.h" #include #include "stir/CPUTimer.h" #include "stir/shared_ptr.h" diff --git a/src/utilities/trim_projdata.cxx b/src/utilities/trim_projdata.cxx index 208f899046..3f6d0d2171 100644 --- a/src/utilities/trim_projdata.cxx +++ b/src/utilities/trim_projdata.cxx @@ -76,8 +76,11 @@ int main(int argc, char **argv) if (in_projdata_ptr->get_proj_data_info_sptr()->get_scanner_ptr()->get_scanner_geometry() == "Cylindrical") { - const ProjDataInfoCylindrical * const in_projdata_info_cyl_ptr = - dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + // const ProjDataInfoCylindrical * const in_projdata_info_cyl_ptr = + // dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + auto proj_data_info_ptr->get_proj_data_info_sptr(); + auto in_projdata_info_cyl_ptr = dynamic_cast(proj_data_info_ptr.get()) + if (in_projdata_info_cyl_ptr== NULL) { error("error converting to cylindrical projection data\n"); @@ -130,8 +133,10 @@ int main(int argc, char **argv) else if (in_projdata_ptr->get_proj_data_info_sptr()->get_scanner_ptr()->get_scanner_geometry() == "BlocksOnCylindrical") { - const ProjDataInfoBlocksOnCylindrical * const in_projdata_info_blk_ptr = - dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + // const ProjDataInfoBlocksOnCylindrical * const in_projdata_info_blk_ptr = + // dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + auto proj_data_info_ptr->get_proj_data_info_sptr(); + auto in_projdata_info_blk_ptr = dynamic_cast(proj_data_info_ptr.get()) if (in_projdata_info_blk_ptr== NULL) { error("error converting to BlocksOnCylindrical projection data\n"); @@ -186,8 +191,10 @@ int main(int argc, char **argv) else if (in_projdata_ptr->get_proj_data_info_sptr()->get_scanner_ptr()->get_scanner_geometry() == "Generic") { - const ProjDataInfoGeneric * const in_projdata_info_gen_ptr = - dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + // const ProjDataInfoGeneric * const in_projdata_info_gen_ptr = + // dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + auto proj_data_info_ptr->get_proj_data_info_sptr(); + auto in_projdata_info_gen_ptr = dynamic_cast(proj_data_info_ptr.get()) if (in_projdata_info_gen_ptr== NULL) { error("error converting to Generic projection data\n"); From ec7a66c63dc8db21da10c9ea5b1bbfe2d449107d Mon Sep 17 00:00:00 2001 From: ssaghamanesh Date: Mon, 9 Sep 2024 12:13:04 +0200 Subject: [PATCH 10/19] Fixed dynamic_cast problem in two utilities --- src/utilities/convert_projdata_types.cxx | 2 +- src/utilities/trim_projdata.cxx | 19 +++++++++++++------ 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/src/utilities/convert_projdata_types.cxx b/src/utilities/convert_projdata_types.cxx index dbbac8f157..d254da003c 100644 --- a/src/utilities/convert_projdata_types.cxx +++ b/src/utilities/convert_projdata_types.cxx @@ -28,7 +28,7 @@ #include "stir/GeometryBlocksOnCylindrical.h" #include "stir/DetectionPosition.h" #include "stir/CartesianCoordinate3D.h" -#include "stir/listmode/DetectorCoordinateMapFromFile.h" +#include "stir/DetectorCoordinateMap.h" #include #include "stir/CPUTimer.h" #include "stir/shared_ptr.h" diff --git a/src/utilities/trim_projdata.cxx b/src/utilities/trim_projdata.cxx index 208f899046..4bac93a7f1 100644 --- a/src/utilities/trim_projdata.cxx +++ b/src/utilities/trim_projdata.cxx @@ -76,8 +76,11 @@ int main(int argc, char **argv) if (in_projdata_ptr->get_proj_data_info_sptr()->get_scanner_ptr()->get_scanner_geometry() == "Cylindrical") { - const ProjDataInfoCylindrical * const in_projdata_info_cyl_ptr = - dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + // const ProjDataInfoCylindrical * const in_projdata_info_cyl_ptr = + // dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + auto proj_data_info_ptr = in_projdata_ptr->get_proj_data_info_sptr(); + auto in_projdata_info_cyl_ptr = dynamic_cast(proj_data_info_ptr.get()); + if (in_projdata_info_cyl_ptr== NULL) { error("error converting to cylindrical projection data\n"); @@ -130,8 +133,10 @@ int main(int argc, char **argv) else if (in_projdata_ptr->get_proj_data_info_sptr()->get_scanner_ptr()->get_scanner_geometry() == "BlocksOnCylindrical") { - const ProjDataInfoBlocksOnCylindrical * const in_projdata_info_blk_ptr = - dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + // const ProjDataInfoBlocksOnCylindrical * const in_projdata_info_blk_ptr = + // dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + auto proj_data_info_ptr = in_projdata_ptr->get_proj_data_info_sptr(); + auto in_projdata_info_blk_ptr = dynamic_cast(proj_data_info_ptr.get()); if (in_projdata_info_blk_ptr== NULL) { error("error converting to BlocksOnCylindrical projection data\n"); @@ -186,8 +191,10 @@ int main(int argc, char **argv) else if (in_projdata_ptr->get_proj_data_info_sptr()->get_scanner_ptr()->get_scanner_geometry() == "Generic") { - const ProjDataInfoGeneric * const in_projdata_info_gen_ptr = - dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + // const ProjDataInfoGeneric * const in_projdata_info_gen_ptr = + // dynamic_cast(in_projdata_ptr->get_proj_data_info_sptr()); + auto proj_data_info_ptr = in_projdata_ptr->get_proj_data_info_sptr(); + auto in_projdata_info_gen_ptr = dynamic_cast(proj_data_info_ptr.get()); if (in_projdata_info_gen_ptr== NULL) { error("error converting to Generic projection data\n"); From a6ad92acdf2360dec4d13fc3bcaed9197ef27fab Mon Sep 17 00:00:00 2001 From: ssaghamanesh Date: Thu, 5 Sep 2024 15:27:32 +0200 Subject: [PATCH 11/19] Initial commit --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 075f775d2a..51c7038fde 100644 --- a/.gitignore +++ b/.gitignore @@ -71,6 +71,7 @@ install_manifest.txt build/ install/ .vscode/ +build_safir_updates/ # MATLAB From 716f95dfe821679af238ba1570e220ce28f636e1 Mon Sep 17 00:00:00 2001 From: ssaghamanesh Date: Thu, 5 Sep 2024 15:28:18 +0200 Subject: [PATCH 12/19] Modified min/max_tang_pos_num to be consistent with ProjDataInfo.cxx --- src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx | 4 ++-- src/buildblock/ProjDataInfoGenericNoArcCorr.cxx | 4 ++-- src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx | 4 +++- 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx index 93b2e49169..a11543cd7a 100644 --- a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx @@ -178,8 +178,8 @@ ProjDataInfoCylindricalNoArcCorr::initialise_uncompressed_view_tangpos_to_det1de assert(fabs(get_phi(Bin(0, 0, 0, 0)) - v_offset) < 1.E-4); assert(fabs(get_phi(Bin(0, get_num_views(), 0, 0)) - v_offset - _PI) < 1.E-4); #endif - const int min_tang_pos_num = -(num_detectors / 2) + 1; - const int max_tang_pos_num = -(num_detectors / 2) + num_detectors; + const int min_tang_pos_num = -(num_detectors / 2); + const int max_tang_pos_num = -(num_detectors / 2) + num_detectors - 1; if (this->get_min_tangential_pos_num() < min_tang_pos_num || this->get_max_tangential_pos_num() > max_tang_pos_num) { diff --git a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx index 6b4ea7de48..930e3be73a 100644 --- a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx @@ -149,8 +149,8 @@ ProjDataInfoGenericNoArcCorr::initialise_uncompressed_view_tangpos_to_det1det2() const int num_detectors = get_scanner_ptr()->get_num_detectors_per_ring(); assert(num_detectors % 2 == 0); - const int min_tang_pos_num = -(num_detectors / 2) + 1; - const int max_tang_pos_num = -(num_detectors / 2) + num_detectors; + const int min_tang_pos_num = -(num_detectors / 2); + const int max_tang_pos_num = -(num_detectors / 2) + num_detectors - 1; if (this->get_min_tangential_pos_num() < min_tang_pos_num || this->get_max_tangential_pos_num() > max_tang_pos_num) { diff --git a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx index 3d596800bb..5613e3fe57 100644 --- a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx +++ b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx @@ -273,9 +273,11 @@ ProjMatrixByBinUsingRayTracing::set_up( voxel_size = image_info_ptr->get_voxel_size(); origin = image_info_ptr->get_origin(); - if (std::abs(origin.x()) > .05F || std::abs(origin.y()) > .05F) + if (std::abs(origin.x()) > .05F || std::abs(origin.y()) > .05F){ error("ProjMatrixByBinUsingRayTracing sadly doesn't support shifted x/y origin yet"); image_info_sptr->get_regular_range(min_index, max_index); + std::cout << "OriginX : " << origin.x() << " ; OriginY : " << origin.y() << "\n"; + } symmetries_sptr.reset(new DataSymmetriesForBins_PET_CartesianGrid(proj_data_info_sptr, density_info_sptr_v, From df7be6a18a510c7f1bd82368261e7fa3786af02c Mon Sep 17 00:00:00 2001 From: ssaghamanesh Date: Thu, 5 Sep 2024 15:27:32 +0200 Subject: [PATCH 13/19] Initial commit --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 075f775d2a..51c7038fde 100644 --- a/.gitignore +++ b/.gitignore @@ -71,6 +71,7 @@ install_manifest.txt build/ install/ .vscode/ +build_safir_updates/ # MATLAB From ac631b5c19fe0178e973f2ad99ddc70e6adebe82 Mon Sep 17 00:00:00 2001 From: ssaghamanesh Date: Thu, 5 Sep 2024 15:28:18 +0200 Subject: [PATCH 14/19] Modified min/max_tang_pos_num to be consistent with ProjDataInfo.cxx --- src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx | 4 ++-- src/buildblock/ProjDataInfoGenericNoArcCorr.cxx | 4 ++-- src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx | 4 +++- 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx index 93b2e49169..a11543cd7a 100644 --- a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx @@ -178,8 +178,8 @@ ProjDataInfoCylindricalNoArcCorr::initialise_uncompressed_view_tangpos_to_det1de assert(fabs(get_phi(Bin(0, 0, 0, 0)) - v_offset) < 1.E-4); assert(fabs(get_phi(Bin(0, get_num_views(), 0, 0)) - v_offset - _PI) < 1.E-4); #endif - const int min_tang_pos_num = -(num_detectors / 2) + 1; - const int max_tang_pos_num = -(num_detectors / 2) + num_detectors; + const int min_tang_pos_num = -(num_detectors / 2); + const int max_tang_pos_num = -(num_detectors / 2) + num_detectors - 1; if (this->get_min_tangential_pos_num() < min_tang_pos_num || this->get_max_tangential_pos_num() > max_tang_pos_num) { diff --git a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx index 6b4ea7de48..930e3be73a 100644 --- a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx @@ -149,8 +149,8 @@ ProjDataInfoGenericNoArcCorr::initialise_uncompressed_view_tangpos_to_det1det2() const int num_detectors = get_scanner_ptr()->get_num_detectors_per_ring(); assert(num_detectors % 2 == 0); - const int min_tang_pos_num = -(num_detectors / 2) + 1; - const int max_tang_pos_num = -(num_detectors / 2) + num_detectors; + const int min_tang_pos_num = -(num_detectors / 2); + const int max_tang_pos_num = -(num_detectors / 2) + num_detectors - 1; if (this->get_min_tangential_pos_num() < min_tang_pos_num || this->get_max_tangential_pos_num() > max_tang_pos_num) { diff --git a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx index 3d596800bb..5613e3fe57 100644 --- a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx +++ b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx @@ -273,9 +273,11 @@ ProjMatrixByBinUsingRayTracing::set_up( voxel_size = image_info_ptr->get_voxel_size(); origin = image_info_ptr->get_origin(); - if (std::abs(origin.x()) > .05F || std::abs(origin.y()) > .05F) + if (std::abs(origin.x()) > .05F || std::abs(origin.y()) > .05F){ error("ProjMatrixByBinUsingRayTracing sadly doesn't support shifted x/y origin yet"); image_info_sptr->get_regular_range(min_index, max_index); + std::cout << "OriginX : " << origin.x() << " ; OriginY : " << origin.y() << "\n"; + } symmetries_sptr.reset(new DataSymmetriesForBins_PET_CartesianGrid(proj_data_info_sptr, density_info_sptr_v, From dacf2943b2b1649dc64684ddd6f82cfe2baf9c80 Mon Sep 17 00:00:00 2001 From: "ssaghamanesh@phys.ethz.ch" Date: Mon, 10 Feb 2025 12:11:56 +0100 Subject: [PATCH 15/19] Added new NCF calculator script. --- .gitignore | 1 + src/utilities/CMakeLists.txt | 1 + ...normfactors_from_cylinder_data_SafirII.cxx | 363 ++++++++++++++++++ 3 files changed, 365 insertions(+) create mode 100644 src/utilities/find_normfactors_from_cylinder_data_SafirII.cxx diff --git a/.gitignore b/.gitignore index 51c7038fde..86c96d4a71 100644 --- a/.gitignore +++ b/.gitignore @@ -72,6 +72,7 @@ build/ install/ .vscode/ build_safir_updates/ +build_tests/ # MATLAB diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt index 1d99a99dfa..b24a7cb2f1 100644 --- a/src/utilities/CMakeLists.txt +++ b/src/utilities/CMakeLists.txt @@ -70,6 +70,7 @@ set(${dir_EXE_SOURCES} find_ML_normfactors.cxx find_ML_singles_from_delayed.cxx find_normfactors_from_cylinder_data.cxx + find_normfactors_from_cylinder_data_SafirII.cxx find_recovery_coefficients_in_image_quality_phantom_nema_nu4.cxx write_sinogram_to_txt.cxx find_sum_projection_of_viewgram_and_sinogram.cxx diff --git a/src/utilities/find_normfactors_from_cylinder_data_SafirII.cxx b/src/utilities/find_normfactors_from_cylinder_data_SafirII.cxx new file mode 100644 index 0000000000..0326833f1a --- /dev/null +++ b/src/utilities/find_normfactors_from_cylinder_data_SafirII.cxx @@ -0,0 +1,363 @@ +/* + Copyright + */ +/*! + + \file + \ingroup utilities + + \brief Find normalisation factors given projection data of a cylinder (component-based method) + + \author Parisa Khateri + +*/ +//#include "stir/Scanner.h" +//#include "stir/stream.h" +//#include "stir/CPUTimer.h" +#include "stir/utilities.h" +#include +#include +#include +#include +#include "stir/shared_ptr.h" +#include "stir/ProjData.h" +#include "stir/CartesianCoordinate3D.h" +#include "stir/LORCoordinates.h" +#include "stir/Bin.h" +#include "stir/ProjDataInterfile.h" +#include "stir/IO/interfile.h" +#include +#include + +USING_NAMESPACE_STIR + + +std::tuple get_crystalTangPos_from_LOR(int view_nr, int tang_pos) { + return {(view_nr + 180 + ((tang_pos - (abs(tang_pos) % 2)) / 2))%180, (view_nr + 270 - ((tang_pos + (abs(tang_pos) % 2)) / 2))%180}; +} + +std::tuple get_crystalAxPos_from_LOR(int segment_nr, int ax_pos) { + if (segment_nr < 0) + { + return {ax_pos - segment_nr, ax_pos}; + } + else + { + return {ax_pos, ax_pos + segment_nr}; + } +} + + +int main(int argc, char **argv) +{ + if (argc!=4) + {//TODO ask the user to give the activity + std::cerr << "Usage: "<< argv[0] + <<" output_file_name_prefix cylider_measured_data cylinder_radius(mm)\n" + <<"only cylinder data are supported. The radius should be the radius of the measured cylinder data.\n" + <<"warning: mind the input order\n"; + return EXIT_FAILURE; + } + + shared_ptr cylinder_projdata_ptr = ProjData::read_from_file(argv[2]); + const std::string output_file_name = argv[1]; + const float R = atof(argv[3]); // cylinder radius + if (R==0) + { + std::cerr << " Radius must be a float value\n" + <<"Usage: "<< argv[0] + <<" output_file_name_prefix cylider_measured_data cylinder_radius\n" + <<"warning: mind the input order\n"; + return EXIT_FAILURE; + } + + //output file + shared_ptr cylinder_pdi_ptr(cylinder_projdata_ptr->get_proj_data_info_ptr()->clone()); + + ProjDataInterfile output_projdata(cylinder_projdata_ptr->get_exam_info_sptr(), cylinder_pdi_ptr, output_file_name); + write_basic_interfile_PDFS_header(output_file_name, output_projdata); + + CartesianCoordinate3D c1, c2; + LORInAxialAndNoArcCorrSinogramCoordinates lor; + + // first find the average number of counts per LOR + std::cout << "Finding the Average NCounts per LOR \n"; + float total_count = 0; + float min_count = std::numeric_limits::max(); // minimum number of counts per LOR + float average_count = 0; //average number of counts per LOR in the active region + int num_active_LORs = 0; //number of LORs which pass through the cylinder + + const int nTang_pos_max = 180; + const int nAx_pos_max = 64; + + for (int seg = cylinder_projdata_ptr->get_min_segment_num(); seg <=cylinder_projdata_ptr->get_max_segment_num(); ++seg) + { + for (int view =cylinder_projdata_ptr->get_min_view_num(); view <=cylinder_projdata_ptr->get_max_view_num(); ++view) + { + Viewgram cylinder_viewgram = cylinder_projdata_ptr->get_viewgram(view, seg); + + for (int ax = cylinder_projdata_ptr->get_min_axial_pos_num(seg); ax <=cylinder_projdata_ptr->get_max_axial_pos_num(seg); ++ax) + { + for (int tang = cylinder_projdata_ptr->get_min_tangential_pos_num(); tang <=cylinder_projdata_ptr->get_max_tangential_pos_num(); ++tang) + { + Bin bin(seg, view, ax, tang); + cylinder_projdata_ptr->get_proj_data_info_ptr()->get_LOR(lor, bin); + LORAs2Points lor_as2points(lor); + LORAs2Points intersection_coords; + if (find_LOR_intersections_with_cylinder(intersection_coords, lor_as2points, R) ==Succeeded::yes) + { //this only succeeds if LOR is intersecting with the cylinder + float N_lor = cylinder_viewgram[ax][tang]; //counts seen by this lor + c1 = intersection_coords.p1(); + c2 = intersection_coords.p2(); + float c12 = sqrt( pow(c1.z()-c2.z(), 2) + pow(c1.y()-c2.y(), 2) + pow(c1.x()-c2.x(), 2) ); // length of intersection of lor with the cylinder + if (c12>0.5) // if LOR intersection is lager than 0.5 mm, check the count per LOR + { + float N_lor_corrected=N_lor/c12; // corrected for the length + total_count+=N_lor_corrected; + num_active_LORs+=1; + if (N_lor_corrected, nAx_pos_max> det_eff {{}}; + std::array, nAx_pos_max> n_alive_channels {{}}; + + std::cout << "Finding the Detector Efficiencies \n"; + for (int view =cylinder_projdata_ptr->get_min_view_num(); view <=cylinder_projdata_ptr->get_max_view_num(); ++view) + { + int seg = 0; + Viewgram cylinder_viewgram = cylinder_projdata_ptr->get_viewgram(view, seg); + for (int ax = cylinder_projdata_ptr->get_min_axial_pos_num(seg); ax <=cylinder_projdata_ptr->get_max_axial_pos_num(seg); ++ax) + { + for (int tang = cylinder_projdata_ptr->get_min_tangential_pos_num(); tang <=cylinder_projdata_ptr->get_max_tangential_pos_num(); ++tang) + { + Bin bin(seg, view, ax, tang); + cylinder_projdata_ptr->get_proj_data_info_ptr()->get_LOR(lor, bin); + float N_lor = cylinder_viewgram[ax][tang]; //counts seen by this lor + if (N_lor > 1) + { + std::tuple crystal_axial = get_crystalAxPos_from_LOR(seg, ax); + std::tuple crystal_tang = get_crystalTangPos_from_LOR(view, tang); + det_eff[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)] += N_lor; + det_eff[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)] += N_lor; + n_alive_channels[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)]++; + n_alive_channels[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)]++; + } + } + } + } + + const int nSectors = 12; + float total_alive_channels = 0; + for (int i = 0; i < (nTang_pos_max/nSectors); i++) { + int dead_channels = 0; + float sum = 0; + for (int j = 0; j < nAx_pos_max; j++) { + for (int k = i; k < nTang_pos_max; k+=(nTang_pos_max/nSectors)) { + sum += det_eff[j][k]; + total_alive_channels += n_alive_channels[j][k]; + } + } + sum = sum / float(nSectors * nAx_pos_max); + total_alive_channels = total_alive_channels / float(nSectors * nAx_pos_max); + for (int j = 0; j < nAx_pos_max; j++) { + for (int k = i; k < nTang_pos_max; k+=(nTang_pos_max/nSectors)) { + float new_val = (det_eff[j][k] * total_alive_channels) / (sum * n_alive_channels[j][k]); + if (det_eff[j][k] < 10.0) { + std::cout << det_eff[j][k] << std::endl; + new_val = 0.00; + } + det_eff[j][k] = new_val; + } + } + } + + float mean_deteff = 0.0; + for (int j = 0; j < nAx_pos_max; j++) { + for (int k = 0; k < nTang_pos_max; k++) { + mean_deteff+=det_eff[j][j]; + } + } + std::cout << "Mean Detector Efficiency" << mean_deteff/float(nAx_pos_max*nTang_pos_max) << std::endl; + + std::cout << "Apply detector efficiencies\n"; + + total_count = 0; + min_count = std::numeric_limits::max(); // minimum number of counts per LOR + average_count = 0; //average number of counts per LOR in the active region + num_active_LORs = 0; //number of LORs which pass through the cylinder + + for (int seg = cylinder_projdata_ptr->get_min_segment_num(); seg <=cylinder_projdata_ptr->get_max_segment_num(); ++seg) + { + for (int view =cylinder_projdata_ptr->get_min_view_num(); view <=cylinder_projdata_ptr->get_max_view_num(); ++view) + { + Viewgram cylinder_viewgram = cylinder_projdata_ptr->get_viewgram(view, seg); + + for (int ax = cylinder_projdata_ptr->get_min_axial_pos_num(seg); ax <=cylinder_projdata_ptr->get_max_axial_pos_num(seg); ++ax) + { + for (int tang = cylinder_projdata_ptr->get_min_tangential_pos_num(); tang <=cylinder_projdata_ptr->get_max_tangential_pos_num(); ++tang) + { + Bin bin(seg, view, ax, tang); + cylinder_projdata_ptr->get_proj_data_info_ptr()->get_LOR(lor, bin); + LORAs2Points lor_as2points(lor); + LORAs2Points intersection_coords; + if (find_LOR_intersections_with_cylinder(intersection_coords, lor_as2points, R) ==Succeeded::yes) + { //this only succeeds if LOR is intersecting with the cylinder + float N_lor = cylinder_viewgram[ax][tang]; //counts seen by this lor + c1 = intersection_coords.p1(); + c2 = intersection_coords.p2(); + float c12 = sqrt( pow(c1.z()-c2.z(), 2) + pow(c1.y()-c2.y(), 2) + pow(c1.x()-c2.x(), 2) ); // length of intersection of lor with the cylinder + if (c12>0.5) // if LOR intersection is lager than 0.5 mm, check the count per LOR + { + std::tuple crystal_axial = get_crystalAxPos_from_LOR(seg, ax); + std::tuple crystal_tang = get_crystalTangPos_from_LOR(view, tang); + float N_lor_corrected=0.0; // corrected for the length and detector efficiencies + if (det_eff[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)] * det_eff[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)] != 0) + { + N_lor_corrected=N_lor/(c12 * det_eff[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)] * det_eff[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)]); // corrected for the length and detector efficiencies + num_active_LORs+=1; + } + total_count+=N_lor_corrected; + if (N_lor_correctedget_max_tangential_pos_num() - cylinder_projdata_ptr->get_min_tangential_pos_num())) + {//TODO ask the user to give the activity + std::cerr << "Error: Geometry Hardcoded for SAFIR-II. Check number of Tangential positions in code to projdata \n"; + return EXIT_FAILURE; + } + if ((nSeg_proj-1)!=cylinder_projdata_ptr->get_max_segment_num()) + {//TODO ask the user to give the activity + std::cerr << "Error: Geometry Hardcoded for SAFIR-II. Check number of segments in code to projdata \n"; + return EXIT_FAILURE; + } + + std::array, nView_proj>, nSeg_proj> geom_coeff {{{}}}; + + for (int seg = 0; seg < nSeg_proj; ++seg) + { + for (int tang = -(nTang_pos_proj/2); tang < (nTang_pos_proj/2); ++tang) + { + for (int base_view = 0; base_view < (nView_proj / (nSectors / 2)); ++base_view) + { + num_active_LORs = 0; + float sum = 0; + for (int true_view = base_view; true_view < nView_proj; true_view+=(nView_proj / (nSectors / 2))) + { + int valid_segments[2] = {seg, -seg}; + for (int true_seg : valid_segments) + { + Viewgram cylinder_viewgram = cylinder_projdata_ptr->get_viewgram(true_view, true_seg); + for (int ax = cylinder_projdata_ptr->get_min_axial_pos_num(true_seg); ax <= cylinder_projdata_ptr->get_max_axial_pos_num(true_seg); ++ax) + { + Bin bin(true_seg, true_view, ax, tang); + cylinder_projdata_ptr->get_proj_data_info_ptr()->get_LOR(lor, bin); + LORAs2Points lor_as2points(lor); + LORAs2Points intersection_coords; + + if (find_LOR_intersections_with_cylinder(intersection_coords, lor_as2points, R) ==Succeeded::yes) + { //this only succeeds if LOR is intersecting with the cylinder + + float N_lor = cylinder_viewgram[ax][tang]; //counts seen by this lor + c1 = intersection_coords.p1(); + c2 = intersection_coords.p2(); + float c12 = sqrt( pow(c1.z()-c2.z(), 2) + pow(c1.y()-c2.y(), 2) + pow(c1.x()-c2.x(), 2) ); // length of intersection of lor with the cylinder + + if (c12>0.5) // if LOR intersection is lager than 0.5 mm, check the count per LOR + { + if (N_lor>1) //if value is large enough + { + std::tuple crystal_axial = get_crystalAxPos_from_LOR(true_seg, ax); + std::tuple crystal_tang = get_crystalTangPos_from_LOR(true_view, tang); + float N_lor_corrected=0.0; // corrected for the length and detector efficiencies + if (det_eff[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)] * det_eff[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)] != 0) + { + float N_lor_corrected=N_lor/(c12 * det_eff[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)] * det_eff[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)]); // corrected for the length and detector efficiencies + num_active_LORs++; + sum += N_lor_corrected; + } + } + } + } + } + } + } + + for (int true_view = base_view; true_view < nView_proj; true_view+=(nView_proj / (nSectors / 2))) + { + if (sum == 0) + { + geom_coeff[seg][true_view][tang+(nTang_pos_proj/2)] = 0.0000; + } + else + { + geom_coeff[seg][true_view][tang+(nTang_pos_proj/2)] = average_count / (sum / num_active_LORs); + } + } + } + } + } + + std::cout << "Find Norm Factors per LOR\n"; + + for (int seg = cylinder_projdata_ptr->get_min_segment_num(); seg <= cylinder_projdata_ptr->get_max_segment_num(); ++seg) + { + for (int view =cylinder_projdata_ptr->get_min_view_num(); view <= cylinder_projdata_ptr->get_max_view_num(); ++view) + { + Viewgram out_viewgram = cylinder_projdata_ptr->get_empty_viewgram(view, seg); + for (int ax = cylinder_projdata_ptr->get_min_axial_pos_num(seg); ax <= cylinder_projdata_ptr->get_max_axial_pos_num(seg); ++ax) + { + for (int tang = cylinder_projdata_ptr->get_min_tangential_pos_num(); tang <= cylinder_projdata_ptr->get_max_tangential_pos_num(); ++tang) + { + Bin bin(seg, view, ax, tang); + cylinder_projdata_ptr->get_proj_data_info_ptr()->get_LOR(lor, bin); + + std::tuple crystal_axial = get_crystalAxPos_from_LOR(seg, ax); + std::tuple crystal_tang = get_crystalTangPos_from_LOR(view, tang); + if (geom_coeff[abs(seg)][view][tang+(nTang_pos_proj/2)] == 0 || det_eff[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)] == 0 || det_eff[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)] == 0) + { + out_viewgram[ax][tang] = 0.0; + } + else + { + out_viewgram[ax][tang] = geom_coeff[abs(seg)][view][tang+(nTang_pos_proj/2)] / (det_eff[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)] * det_eff[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)]); + } + } + } + output_projdata.set_viewgram(out_viewgram); + } + } +} From 8959547061910dc603b0963bec7c82e869714ce2 Mon Sep 17 00:00:00 2001 From: "ssaghamanesh@phys.ethz.ch" Date: Mon, 10 Feb 2025 12:13:22 +0100 Subject: [PATCH 16/19] .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 51c7038fde..86c96d4a71 100644 --- a/.gitignore +++ b/.gitignore @@ -72,6 +72,7 @@ build/ install/ .vscode/ build_safir_updates/ +build_tests/ # MATLAB From c6e847a4763a6336c2a45d6dd7cf18fa26aaf67d Mon Sep 17 00:00:00 2001 From: "ssaghamanesh@phys.ethz.ch" Date: Mon, 10 Feb 2025 12:57:12 +0100 Subject: [PATCH 17/19] .gitignore --- .gitignore | 2 -- 1 file changed, 2 deletions(-) diff --git a/.gitignore b/.gitignore index 86c96d4a71..075f775d2a 100644 --- a/.gitignore +++ b/.gitignore @@ -71,8 +71,6 @@ install_manifest.txt build/ install/ .vscode/ -build_safir_updates/ -build_tests/ # MATLAB From b5f4f78057ae9546dc0a96a11faf893d97edfe53 Mon Sep 17 00:00:00 2001 From: "ssaghamanesh@phys.ethz.ch" Date: Mon, 10 Feb 2025 13:04:08 +0100 Subject: [PATCH 18/19] Removed a typo. --- src/buildblock/Scanner.cxx | 58 -------------------------------------- 1 file changed, 58 deletions(-) diff --git a/src/buildblock/Scanner.cxx b/src/buildblock/Scanner.cxx index dd42a1f6f7..fbceb63906 100644 --- a/src/buildblock/Scanner.cxx +++ b/src/buildblock/Scanner.cxx @@ -1146,64 +1146,6 @@ Scanner::Scanner(Type scanner_type) ); break; - case SafirI: - set_params(SafirI, string_list("SafirI"), - 24, //num_rings_v - 150, //max_num_non_arccorrected_bins_v, - 150, //default_num_arccorrected_bins_v, - 180, //num_detectors_per_ring_v - 64.05, // inner_ring_radius_v - 5, //average_depth_of_interaction_v - 2.2, //ring_spacing_v - 1.1, //bin_size_v - 0, //intrinsic_tilt_v - 3, //num_axial_blocks_per_bucket_v - 1, //num_transaxial_blocks_per_bucket_v - 8, //num_axial_crystals_per_block_v - 15, //num_transaxial_crystals_per_block_v - 1, //num_axial_crystals_per_singles_unit_v - 1, //num_transaxial_crystals_per_singles_unit_v - 1, //num_detector_layers_v - 0.11, //energy_resolution_v - 511, //reference_energy_v - "", //scanner_geometry_v - 2.2, //axial_crystal_spacing_v - 2.2, //transaxial_crystal_spacing_v - 18.1, //axial_block_spacing_v - 33.6, //transaxial_block_spacing_v - ""//crystal_map_file_name_v - ); - break; - - case SafirII: - set_params(SafirII, string_list("SafirII"), - 64, //num_rings_v - 150, //max_num_non_arccorrected_bins_v, - 150, //default_num_arccorrected_bins_v, - 180, //num_detectors_per_ring_v - 64.05, // inner_ring_radius_v - 5, //average_depth_of_interaction_v - 2.2, //ring_spacing_v - 1.1, //bin_size_v - 0, //intrinsic_tilt_v - 8, //num_axial_blocks_per_bucket_v - 1, //num_transaxial_blocks_per_bucket_v - 8, //num_axial_crystals_per_block_v - 15, //num_transaxial_crystals_per_block_v - 1, //num_axial_crystals_per_singles_unit_v - 1, //num_transaxial_crystals_per_singles_unit_v - 1, //num_detector_layers_v - 0.11, //energy_resolution_v - 511, //reference_energy_v - "", //scanner_geometry_v - 2.2, //axial_crystal_spacing_v - 2.2, //transaxial_crystal_spacing_v - 18.1, //axial_block_spacing_v - 33.6, //transaxial_block_spacing_v - ""//crystal_map_file_name_v - ); - break; - case UPENN_5rings: set_params(UPENN_5rings, string_list("UPENN_5rings"), From 197a288b42263e4b9aeea442c69b26de4803228b Mon Sep 17 00:00:00 2001 From: "ssaghamanesh@phys.ethz.ch" Date: Mon, 10 Feb 2025 15:35:30 +0100 Subject: [PATCH 19/19] Modified according to stir6. --- .../find_normfactors_from_cylinder_data_SafirII.cxx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/utilities/find_normfactors_from_cylinder_data_SafirII.cxx b/src/utilities/find_normfactors_from_cylinder_data_SafirII.cxx index 0326833f1a..43ef676cde 100644 --- a/src/utilities/find_normfactors_from_cylinder_data_SafirII.cxx +++ b/src/utilities/find_normfactors_from_cylinder_data_SafirII.cxx @@ -72,7 +72,7 @@ int main(int argc, char **argv) } //output file - shared_ptr cylinder_pdi_ptr(cylinder_projdata_ptr->get_proj_data_info_ptr()->clone()); + shared_ptr cylinder_pdi_ptr(cylinder_projdata_ptr->get_proj_data_info_sptr()->clone()); ProjDataInterfile output_projdata(cylinder_projdata_ptr->get_exam_info_sptr(), cylinder_pdi_ptr, output_file_name); write_basic_interfile_PDFS_header(output_file_name, output_projdata); @@ -101,7 +101,7 @@ int main(int argc, char **argv) for (int tang = cylinder_projdata_ptr->get_min_tangential_pos_num(); tang <=cylinder_projdata_ptr->get_max_tangential_pos_num(); ++tang) { Bin bin(seg, view, ax, tang); - cylinder_projdata_ptr->get_proj_data_info_ptr()->get_LOR(lor, bin); + cylinder_projdata_ptr->get_proj_data_info_sptr()->get_LOR(lor, bin); LORAs2Points lor_as2points(lor); LORAs2Points intersection_coords; if (find_LOR_intersections_with_cylinder(intersection_coords, lor_as2points, R) ==Succeeded::yes) @@ -144,7 +144,7 @@ int main(int argc, char **argv) for (int tang = cylinder_projdata_ptr->get_min_tangential_pos_num(); tang <=cylinder_projdata_ptr->get_max_tangential_pos_num(); ++tang) { Bin bin(seg, view, ax, tang); - cylinder_projdata_ptr->get_proj_data_info_ptr()->get_LOR(lor, bin); + cylinder_projdata_ptr->get_proj_data_info_sptr()->get_LOR(lor, bin); float N_lor = cylinder_viewgram[ax][tang]; //counts seen by this lor if (N_lor > 1) { @@ -210,7 +210,7 @@ int main(int argc, char **argv) for (int tang = cylinder_projdata_ptr->get_min_tangential_pos_num(); tang <=cylinder_projdata_ptr->get_max_tangential_pos_num(); ++tang) { Bin bin(seg, view, ax, tang); - cylinder_projdata_ptr->get_proj_data_info_ptr()->get_LOR(lor, bin); + cylinder_projdata_ptr->get_proj_data_info_sptr()->get_LOR(lor, bin); LORAs2Points lor_as2points(lor); LORAs2Points intersection_coords; if (find_LOR_intersections_with_cylinder(intersection_coords, lor_as2points, R) ==Succeeded::yes) @@ -284,7 +284,7 @@ int main(int argc, char **argv) for (int ax = cylinder_projdata_ptr->get_min_axial_pos_num(true_seg); ax <= cylinder_projdata_ptr->get_max_axial_pos_num(true_seg); ++ax) { Bin bin(true_seg, true_view, ax, tang); - cylinder_projdata_ptr->get_proj_data_info_ptr()->get_LOR(lor, bin); + cylinder_projdata_ptr->get_proj_data_info_sptr()->get_LOR(lor, bin); LORAs2Points lor_as2points(lor); LORAs2Points intersection_coords; @@ -343,7 +343,7 @@ int main(int argc, char **argv) for (int tang = cylinder_projdata_ptr->get_min_tangential_pos_num(); tang <= cylinder_projdata_ptr->get_max_tangential_pos_num(); ++tang) { Bin bin(seg, view, ax, tang); - cylinder_projdata_ptr->get_proj_data_info_ptr()->get_LOR(lor, bin); + cylinder_projdata_ptr->get_proj_data_info_sptr()->get_LOR(lor, bin); std::tuple crystal_axial = get_crystalAxPos_from_LOR(seg, ax); std::tuple crystal_tang = get_crystalTangPos_from_LOR(view, tang);