diff --git a/src/buildblock/interpolate_projdata.cxx b/src/buildblock/interpolate_projdata.cxx index b0e36f8e0c..9c7d5778d8 100644 --- a/src/buildblock/interpolate_projdata.cxx +++ b/src/buildblock/interpolate_projdata.cxx @@ -249,16 +249,20 @@ interpolate_projdata(ProjData& proj_data_out, }; } else - { // for BlocksOnCylindrical, views and tangential positions are not subsampled and can be mapped 1:1 - if (proj_data_in_info.get_num_tangential_poss() != proj_data_out_info.get_num_tangential_poss()) + { // for BlocksOnCylindrical, views and tangential positions are scaled by a fixed value that should be a factor + // of the number of detectors per bucket, such that the intersections between buckets are handled correctly + auto scale_factor = (double)proj_data_out_info.get_num_tangential_poss() / (double)proj_data_in_info.get_num_tangential_poss(); + + if (proj_data_in_info.get_scanner_sptr()->get_num_detectors_per_ring() * (int)scale_factor != + proj_data_out_info.get_scanner_sptr()->get_num_detectors_per_ring()) { - error("Interpolation of BlocksOnCylindrical scanners assumes that number of tangential positions " - "is the same in the downsampled scanner."); + warning("The scanner downsampling is not done by an integer factor. This can produce issues in the interpolation due to the block geometry."); } - if (proj_data_in_info.get_num_views() != proj_data_out_info.get_num_views()) + if (proj_data_in_info.get_scanner_sptr()->get_num_transaxial_crystals_per_bucket() * (int)scale_factor != + proj_data_out_info.get_scanner_sptr()->get_num_transaxial_crystals_per_bucket()) { - error("Interpolation of BlocksOnCylindrical scanners assumes that number of views " - "is the same in the downsampled scanner."); + warning("The scanner downsampling per bucket is not done by an integer factor. " + "This can produce issues in the interpolation due to the block geometry"); } // only extending in axial direction - an extension of 2 was found to be sufficient @@ -275,7 +279,7 @@ interpolate_projdata(ProjData& proj_data_out, } // define a function to translate indices in the output proj data to indices in input proj data - index_converter = [&proj_data_out_info, m_offset, m_sampling](const BasicCoordinate<3, int>& index_out) + index_converter = [&proj_data_out_info, m_offset, m_sampling, scale_factor](const BasicCoordinate<3, int>& index_out) -> BasicCoordinate<3, double> { // translate index on output to coordinate @@ -285,8 +289,8 @@ interpolate_projdata(ProjData& proj_data_out, // translate to indices in input proj data BasicCoordinate<3, double> index_in; index_in[1] = (out_m - m_offset) / m_sampling; - index_in[2] = index_out[2]; - index_in[3] = index_out[3]; + index_in[2] = index_out[2] / scale_factor; + index_in[3] = index_out[3] / scale_factor; return index_in; }; diff --git a/src/include/stir/scatter/ScatterEstimation.h b/src/include/stir/scatter/ScatterEstimation.h index 75980d1220..bc4e3ea25b 100644 --- a/src/include/stir/scatter/ScatterEstimation.h +++ b/src/include/stir/scatter/ScatterEstimation.h @@ -172,6 +172,8 @@ class ScatterEstimation: public ParsingObject inline void set_num_iterations(int); + inline void set_downsample_scanner(bool downsample_scanner, int downsampled_detectors_per_ring = 0); + void set_output_scatter_estimate_prefix(const std::string&); void set_export_scatter_estimates_of_each_iteration(bool); @@ -385,6 +387,8 @@ class ScatterEstimation: public ParsingObject float min_scale_value; bool downsample_scanner_bool; + int downsampled_detectors_per_ring; + //! unsigned int half_filter_width; diff --git a/src/include/stir/scatter/ScatterEstimation.inl b/src/include/stir/scatter/ScatterEstimation.inl index ef1260bcf5..5a2f6dbfc0 100644 --- a/src/include/stir/scatter/ScatterEstimation.inl +++ b/src/include/stir/scatter/ScatterEstimation.inl @@ -73,4 +73,12 @@ set_num_iterations(int arg) this->num_scatter_iterations = arg; } +void +ScatterEstimation:: +set_downsample_scanner(bool downsample_scanner, int downsampled_detectors_per_ring) +{ + this->downsample_scanner_bool = downsample_scanner; + this->downsampled_detectors_per_ring = downsampled_detectors_per_ring; +} + END_NAMESPACE_STIR diff --git a/src/scatter_buildblock/ScatterEstimation.cxx b/src/scatter_buildblock/ScatterEstimation.cxx index 50267c9ab9..f1ffa39bcf 100644 --- a/src/scatter_buildblock/ScatterEstimation.cxx +++ b/src/scatter_buildblock/ScatterEstimation.cxx @@ -78,6 +78,7 @@ set_defaults() this->override_scanner_template = true; this->override_density_image = true; this->downsample_scanner_bool = true; + this->downsampled_detectors_per_ring = 0; this->remove_interleaving = true; this->atten_image_filename = ""; this->atten_coeff_filename = ""; @@ -146,6 +147,8 @@ initialise_keymap() &this->scatter_sim_par_filename); this->parser.add_key("use scanner downsampling in scatter simulation", &this->downsample_scanner_bool); + this->parser.add_key("override number of downsampled detectors per ring", + &this->downsampled_detectors_per_ring); this->parser.add_key("override attenuation image", &this->override_density_image); @@ -605,7 +608,7 @@ set_up() } if (this->downsample_scanner_bool) - this->scatter_simulation_sptr->downsample_scanner(); + this->scatter_simulation_sptr->downsample_scanner(-1, this->downsampled_detectors_per_ring); // Check if Load a mask proj_data diff --git a/src/scatter_buildblock/ScatterSimulation.cxx b/src/scatter_buildblock/ScatterSimulation.cxx index 7d28caceda..04b2ad91e6 100644 --- a/src/scatter_buildblock/ScatterSimulation.cxx +++ b/src/scatter_buildblock/ScatterSimulation.cxx @@ -909,39 +909,40 @@ ScatterSimulation::downsample_scanner(int new_num_rings, int new_num_dets) float approx_num_non_arccorrected_bins; if (new_scanner_sptr->get_scanner_geometry()!="Cylindrical") { + if (new_num_dets <= 0) + { // by default, do not downsample the detectors per ring for BlocksOnCylindrical new_num_dets = this->proj_data_info_sptr->get_scanner_ptr()->get_num_detectors_per_ring(); - approx_num_non_arccorrected_bins = this->proj_data_info_sptr->get_num_tangential_poss(); - // preserve the length of the scanner the following includes gaps - float scanner_length_block = new_scanner_sptr->get_num_axial_buckets()* - new_scanner_sptr->get_num_axial_blocks_per_bucket()* - new_scanner_sptr->get_axial_block_spacing(); - new_scanner_sptr->set_num_axial_blocks_per_bucket(1); -// new_scanner_sptr->set_num_transaxial_blocks_per_bucket(1); - - - new_scanner_sptr->set_num_rings(new_num_rings); -// float transaxial_bucket_spacing=old_scanner_ptr->get_transaxial_block_spacing() -// *old_scanner_ptr->get_num_transaxial_blocks_per_bucket(); - float new_ring_spacing=scanner_length_block/new_scanner_sptr->get_num_rings(); -// int num_trans_buckets=old_scanner_ptr->get_num_transaxial_buckets(); -// get a new number of detectors that is a multiple of the number of buckets to preserve scanner shape -// float frac,whole; -// frac = std::modf(float(new_num_dets/new_scanner_sptr->get_num_transaxial_buckets()), &whole); -// int newest_num_dets=whole*new_scanner_sptr->get_num_transaxial_buckets(); -// new_scanner_sptr->set_num_detectors_per_ring(newest_num_dets); -// int new_transaxial_dets_per_bucket=newest_num_dets/num_trans_buckets; -// float new_det_spacing=transaxial_bucket_spacing/new_transaxial_dets_per_bucket; - - new_scanner_sptr->set_axial_crystal_spacing(new_ring_spacing); - new_scanner_sptr->set_ring_spacing(new_ring_spacing); - new_scanner_sptr->set_num_axial_crystals_per_block(new_num_rings); - new_scanner_sptr->set_axial_block_spacing(new_ring_spacing - * new_scanner_sptr->get_num_axial_crystals_per_block()); + } + + // extend the bins by a small amount to avoid edge-effects, at the expense of longer computation time + approx_num_non_arccorrected_bins = ceil(this->proj_data_info_sptr->get_num_tangential_poss() * + float(new_num_dets) / old_scanner_ptr->get_num_detectors_per_ring()) + 1; + // preserve the length of the scanner the following includes gaps + float scanner_length_block = new_scanner_sptr->get_num_axial_buckets()* + new_scanner_sptr->get_num_axial_blocks_per_bucket()* + new_scanner_sptr->get_axial_block_spacing(); + new_scanner_sptr->set_num_axial_blocks_per_bucket(1); + new_scanner_sptr->set_num_transaxial_blocks_per_bucket(1); + + + new_scanner_sptr->set_num_rings(new_num_rings); + float transaxial_bucket_spacing = old_scanner_ptr->get_transaxial_block_spacing() * + old_scanner_ptr->get_num_transaxial_blocks_per_bucket(); + float new_ring_spacing=scanner_length_block/new_scanner_sptr->get_num_rings(); + int num_trans_buckets=old_scanner_ptr->get_num_transaxial_buckets(); + // get a new number of detectors that is a multiple of the number of buckets to preserve scanner shape + new_scanner_sptr->set_num_detectors_per_ring(new_num_dets); + int new_transaxial_dets_per_bucket = new_num_dets / num_trans_buckets; + float new_det_spacing = transaxial_bucket_spacing / new_transaxial_dets_per_bucket; + + new_scanner_sptr->set_axial_crystal_spacing(new_ring_spacing); + new_scanner_sptr->set_ring_spacing(new_ring_spacing); + new_scanner_sptr->set_num_axial_crystals_per_block(new_num_rings); + new_scanner_sptr->set_axial_block_spacing(new_ring_spacing * new_num_rings); -// new_scanner_sptr->set_num_transaxial_crystals_per_block(new_transaxial_dets_per_bucket); -// new_scanner_sptr->set_transaxial_crystal_spacing(new_det_spacing); -// new_scanner_sptr->set_transaxial_block_spacing(new_det_spacing -// * new_scanner_sptr->get_num_transaxial_crystals_per_block()); + new_scanner_sptr->set_num_transaxial_crystals_per_block(new_transaxial_dets_per_bucket); + new_scanner_sptr->set_transaxial_crystal_spacing(new_det_spacing); + new_scanner_sptr->set_transaxial_block_spacing(transaxial_bucket_spacing); } else{ if (new_num_dets <= 0)