From 69486f5ac29c41a9ece446f6b0efda7ee52967ed Mon Sep 17 00:00:00 2001 From: ssaghamanesh Date: Thu, 5 Sep 2024 11:21:07 +0200 Subject: [PATCH] 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; +}