From ac8462c2ba4a45302bb50812c99296c07964eb6e Mon Sep 17 00:00:00 2001 From: "ssaghamanesh@phys.ethz.ch" Date: Mon, 10 Feb 2025 16:24:04 +0100 Subject: [PATCH 1/4] Initial commit: added SafirI/SafirII scanners to Scanner class. --- src/buildblock/Scanner.cxx | 64 ++++++++++++++++++++++++++++++++++++++ src/include/stir/Scanner.h | 2 ++ 2 files changed, 66 insertions(+) diff --git a/src/buildblock/Scanner.cxx b/src/buildblock/Scanner.cxx index 6fbfe6ad9d..0da145d7ca 100644 --- a/src/buildblock/Scanner.cxx +++ b/src/buildblock/Scanner.cxx @@ -1146,6 +1146,70 @@ 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.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 5b2210c694..94759c66ca 100644 --- a/src/include/stir/Scanner.h +++ b/src/include/stir/Scanner.h @@ -175,6 +175,8 @@ class Scanner Allegro, GeminiTF, SAFIRDualRingPrototype, + SafirI, + SafirII, UPENN_5rings, UPENN_5rings_no_gaps, UPENN_6rings, From 8680f5276f375818fbd6e0f3ce401bf336a7822f Mon Sep 17 00:00:00 2001 From: "ssaghamanesh@phys.ethz.ch" Date: Mon, 10 Feb 2025 16:24:57 +0100 Subject: [PATCH 2/4] Adding two utilities for conversion and trimming of sinograms in generic geometry. --- src/utilities/CMakeLists.txt | 65 ++++++ src/utilities/convert_projdata_types.cxx | 95 +++++++++ src/utilities/trim_projdata.cxx | 253 +++++++++++++++++++++++ 3 files changed, 413 insertions(+) create mode 100644 src/utilities/convert_projdata_types.cxx create mode 100644 src/utilities/trim_projdata.cxx diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt index 05faf85500..09423fbc9d 100644 --- a/src/utilities/CMakeLists.txt +++ b/src/utilities/CMakeLists.txt @@ -77,9 +77,74 @@ if (NOT MINI_STIR) stir_timings.cxx ) +<<<<<<< HEAD if (HAVE_ITK) list(APPEND ${dir_EXE_SOURCES} SPECT_dicom_to_interfile.cxx) endif() +======= + manip_image.cxx + stir_write_pgm.cxx + manip_projdata.cxx + display_projdata.cxx + do_linear_regression.cxx + postfilter.cxx + compare_projdata.cxx + compare_image.cxx + extract_segments.cxx + extract_single_images_from_dynamic_image.cxx + estimate_triple_energy_window_scatter_sinogram.cxx + correct_projdata.cxx + stir_math.cxx + list_projdata_info.cxx + list_image_info.cxx + list_image_values.cxx + list_detector_and_bin_info.cxx + find_maxima_in_image.cxx + create_multi_header.cxx + create_projdata_template.cxx + SSRB.cxx + invert_axis.cxx + poisson_noise.cxx + get_time_frame_info.cxx + generate_image.cxx + list_ROI_values.cxx + zoom_image.cxx + find_fwhm_in_image.cxx + abs_image.cxx + convert_to_binary_image.cxx + rebin_projdata.cxx + write_proj_matrix_by_bin.cxx + forward_project.cxx + back_project.cxx + calculate_attenuation_coefficients.cxx + compute_sqrt_Hessian_row_sum.cxx + attenuation_coefficients_to_projections.cxx + conv_GATE_raw_ECAT_projdata_to_interfile.cxx + conv_gipl_to_interfile.cxx + conv_interfile_to_gipl.cxx + shift_image.cxx + stir_config.cxx + stir_list_registries.cxx + shift_image_origin.cxx + warp_and_accumulate_gated_images.cxx + warp_image.cxx + zeropad_planes.cxx + apply_normfactors3D.cxx + apply_normfactors.cxx + construct_randoms_from_singles.cxx + find_ML_normfactors3D.cxx + find_ML_normfactors.cxx + find_ML_singles_from_delayed.cxx + find_normfactors_from_cylinder_data.cxx + find_recovery_coefficients_in_image_quality_phantom_nema_nu4.cxx + write_sinogram_to_txt.cxx + 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 + +>>>>>>> 30780c318 (Adding two utilities for conversion and trimming of sinograms in generic geometry.) if (HAVE_HDF5) list(APPEND ${dir_EXE_SOURCES} diff --git a/src/utilities/convert_projdata_types.cxx b/src/utilities/convert_projdata_types.cxx new file mode 100644 index 0000000000..d254da003c --- /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/DetectorCoordinateMap.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..4bac93a7f1 --- /dev/null +++ b/src/utilities/trim_projdata.cxx @@ -0,0 +1,253 @@ +// +// +/* + +*/ +/*! + + \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()); + 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"); + } + 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()); + 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"); + } + 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()); + 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"); + } + 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 404647c1227fb85fd3d80a6496dde0ba177026e9 Mon Sep 17 00:00:00 2001 From: Somayeh Date: Wed, 12 Nov 2025 12:45:52 +0100 Subject: [PATCH 3/4] Added two utilities to utilities/CMakeLists and edited release_6.4. --- documentation/release_6.4.htm | 6 ++++ src/utilities/CMakeLists.txt | 67 ++--------------------------------- 2 files changed, 8 insertions(+), 65 deletions(-) diff --git a/documentation/release_6.4.htm b/documentation/release_6.4.htm index efdddc2c09..6b918f9d43 100644 --- a/documentation/release_6.4.htm +++ b/documentation/release_6.4.htm @@ -40,6 +40,9 @@

Summary for end users (also to be read by developers)

New functionality

General

    +
  • + Two SAFIR scanners (Safir-I and Safir-II) with a Generic geometry were added to the list of scanners + (see PR for details).

Python

    @@ -47,6 +50,9 @@

    Python

    Utilities

      +
    • + Two utilities were added to convert and trim sinograms in a generic geometry + (see PR for details).

    Changed functionality

    diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt index 09423fbc9d..ce75195fec 100644 --- a/src/utilities/CMakeLists.txt +++ b/src/utilities/CMakeLists.txt @@ -75,76 +75,13 @@ if (NOT MINI_STIR) 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 ) -<<<<<<< HEAD if (HAVE_ITK) list(APPEND ${dir_EXE_SOURCES} SPECT_dicom_to_interfile.cxx) endif() -======= - manip_image.cxx - stir_write_pgm.cxx - manip_projdata.cxx - display_projdata.cxx - do_linear_regression.cxx - postfilter.cxx - compare_projdata.cxx - compare_image.cxx - extract_segments.cxx - extract_single_images_from_dynamic_image.cxx - estimate_triple_energy_window_scatter_sinogram.cxx - correct_projdata.cxx - stir_math.cxx - list_projdata_info.cxx - list_image_info.cxx - list_image_values.cxx - list_detector_and_bin_info.cxx - find_maxima_in_image.cxx - create_multi_header.cxx - create_projdata_template.cxx - SSRB.cxx - invert_axis.cxx - poisson_noise.cxx - get_time_frame_info.cxx - generate_image.cxx - list_ROI_values.cxx - zoom_image.cxx - find_fwhm_in_image.cxx - abs_image.cxx - convert_to_binary_image.cxx - rebin_projdata.cxx - write_proj_matrix_by_bin.cxx - forward_project.cxx - back_project.cxx - calculate_attenuation_coefficients.cxx - compute_sqrt_Hessian_row_sum.cxx - attenuation_coefficients_to_projections.cxx - conv_GATE_raw_ECAT_projdata_to_interfile.cxx - conv_gipl_to_interfile.cxx - conv_interfile_to_gipl.cxx - shift_image.cxx - stir_config.cxx - stir_list_registries.cxx - shift_image_origin.cxx - warp_and_accumulate_gated_images.cxx - warp_image.cxx - zeropad_planes.cxx - apply_normfactors3D.cxx - apply_normfactors.cxx - construct_randoms_from_singles.cxx - find_ML_normfactors3D.cxx - find_ML_normfactors.cxx - find_ML_singles_from_delayed.cxx - find_normfactors_from_cylinder_data.cxx - find_recovery_coefficients_in_image_quality_phantom_nema_nu4.cxx - write_sinogram_to_txt.cxx - 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 - ->>>>>>> 30780c318 (Adding two utilities for conversion and trimming of sinograms in generic geometry.) if (HAVE_HDF5) list(APPEND ${dir_EXE_SOURCES} From 3e946bbde76ebb0ccdde8f1200432eb231a94aa0 Mon Sep 17 00:00:00 2001 From: Somayeh Date: Sun, 28 Jun 2026 17:45:20 +0200 Subject: [PATCH 4/4] Added normalization script for SafirII. --- src/utilities/CMakeLists.txt | 1 + ...normfactors_from_cylinder_data_SafirII.cxx | 363 ++++++++++++++++++ 2 files changed, 364 insertions(+) create mode 100755 src/utilities/find_normfactors_from_cylinder_data_SafirII.cxx diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt index ce75195fec..960697f1e6 100644 --- a/src/utilities/CMakeLists.txt +++ b/src/utilities/CMakeLists.txt @@ -70,6 +70,7 @@ if (NOT MINI_STIR) 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 100755 index 0000000000..43ef676cde --- /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_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); + + 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_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) + { //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_sptr()->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_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) + { //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_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) + { //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_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); + 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); + } + } +}