From 4407e3b396e0d7e5a0fe345d89e06ef20048a5f4 Mon Sep 17 00:00:00 2001 From: ssaghamanesh Date: Thu, 5 Sep 2024 11:20:26 +0200 Subject: [PATCH 1/4] 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 2f48f1438..9ddcab555 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 e5b7332bd..97ee665a3 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 2/4] 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 000000000..dbbac8f15 --- /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 000000000..208f89904 --- /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 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 3/4] 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 42e09b831..1d99a99df 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 ec7a66c63dc8db21da10c9ea5b1bbfe2d449107d Mon Sep 17 00:00:00 2001 From: ssaghamanesh Date: Mon, 9 Sep 2024 12:13:04 +0200 Subject: [PATCH 4/4] 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 dbbac8f15..d254da003 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 208f89904..4bac93a7f 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");