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/buildblock/Scanner.cxx b/src/buildblock/Scanner.cxx index 4cf34c5c6d..fbceb63906 100644 --- a/src/buildblock/Scanner.cxx +++ b/src/buildblock/Scanner.cxx @@ -1120,7 +1120,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 @@ -1134,15 +1134,15 @@ 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; 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, 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, diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt index 42e09b8318..b24a7cb2f1 100644 --- a/src/utilities/CMakeLists.txt +++ b/src/utilities/CMakeLists.txt @@ -70,11 +70,14 @@ 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 separate_true_from_random_scatter_for_necr.cxx stir_timings.cxx + convert_projdata_types.cxx + trim_projdata.cxx ) if (HAVE_ITK) 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/find_normfactors_from_cylinder_data_SafirII.cxx b/src/utilities/find_normfactors_from_cylinder_data_SafirII.cxx new file mode 100644 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); + } + } +} 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; +}