diff --git a/recon_test_pack/DFM3D_test_sim.par b/recon_test_pack/DFM3D_test_sim.par new file mode 100644 index 000000000..1ffd86b9d --- /dev/null +++ b/recon_test_pack/DFM3D_test_sim.par @@ -0,0 +1,13 @@ +DFM3DParameters := +; test file for DFM3D +input file := my_precorrected_sino${suffix}.hs + +xy output image size (in pixels) := -1 +zoom := 0.5 + +noise filter := 0 + +output filename prefix := my_test_sim_image_DFM3D + + +end:= diff --git a/recon_test_pack/run_test_simulate_and_recon.sh b/recon_test_pack/run_test_simulate_and_recon.sh index 2eb7a1664..33c5a85e6 100755 --- a/recon_test_pack/run_test_simulate_and_recon.sh +++ b/recon_test_pack/run_test_simulate_and_recon.sh @@ -61,6 +61,7 @@ echo "Using `command -v SRT2D`" echo "Using `command -v SRT2DSPECT`" echo "Using `command -v GRD2D`" echo "Using `command -v DDSR2D`" +echo "Using `command -v DFM3D`" # first need to set this to the C locale, as this is what the STIR utilities use # otherwise, awk might interpret floating point numbers incorrectly @@ -97,6 +98,14 @@ if [ $? -ne 0 ]; then exit 1 fi +## 3D data (for DFM3D test) +D3_suffix=_D3 +./simulate_data_for_tests.sh --D3 --suffix "$D3_suffix" +if [ $? -ne 0 ]; then + echo "Error running simulation" + exit 1 +fi + error_log_files="" input_image=my_uniform_cylinder.hv @@ -109,7 +118,7 @@ input_ROI_mean=`awk 'NR>2 {print $2}' ${input_image}.roistats` # warning: currently OSMAPOSL needs to be run before OSSPS as # the OSSPS par file uses an OSMAPOSL result as initial image # and reuses its subset sensitivities -for recon in FBP2D FBP3DRP SRT2D SRT2DSPECT GRD2D DDSR2D OSMAPOSL OSSPS ; do +for recon in FBP2D FBP3DRP SRT2D SRT2DSPECT GRD2D DDSR2D DFM3D OSMAPOSL OSSPS ; do echo "========== Testing `command -v ${recon}`" # Check if we have CUDA code and parallelproj. # If so, check for test files in CUDA/* @@ -140,10 +149,12 @@ for recon in FBP2D FBP3DRP SRT2D SRT2DSPECT GRD2D DDSR2D OSMAPOSL OSSPS ; do is_analytic=1 elif expr "$recon" : DDSR > /dev/null; then is_analytic=1 + elif expr "$recon" : DFM > /dev/null; then + is_analytic=1 fi if [ $is_analytic = 1 ]; then if expr "$dataSuffix" : '.*TOF.*' > /dev/null; then - echo "Skipping TOF as not yet supported for FBP, SRT, GRD and DDSR." + echo "Skipping TOF as not yet supported for FBP, SRT, GRD, DDSR and DFM." break fi if expr "$recon" : SRT2DSPECT > /dev/null; then @@ -152,6 +163,16 @@ for recon in FBP2D FBP3DRP SRT2D SRT2DSPECT GRD2D DDSR2D OSMAPOSL OSSPS ; do elif expr "$recon" : DDSR2D > /dev/null; then suffix=$SPECT_suffix export suffix + elif expr "$recon" : DFM3D > /dev/null; then + suffix=$D3_suffix + export suffix + echo "Running precorrection" + correct_projdata correct_projdata_simulation.par > my_correct_projdata_simulation.log 2>&1 + if [ $? -ne 0 ]; then + echo "Error running precorrection. CHECK my_correct_projdata_simulation.log" + error_log_files="${error_log_files} my_correct_projdata_simulation.log" + break + fi else suffix=$zero_view_suffix export suffix diff --git a/recon_test_pack/simulate_data_for_tests.sh b/recon_test_pack/simulate_data_for_tests.sh index 18e04f42e..a08bb6de1 100755 --- a/recon_test_pack/simulate_data_for_tests.sh +++ b/recon_test_pack/simulate_data_for_tests.sh @@ -30,6 +30,7 @@ force_zero_view_offset=0 TOF=0 suffix="" SPECT=0 +D3=0 # # Parse option arguments (--) # Note that the -- is required to suppress interpretation of $1 as options @@ -47,6 +48,9 @@ do elif test "$1" = "--SPECT" then SPECT=1 + elif test "$1" = "--D3" + then + D3=1 elif test "$1" = "--suffix" then suffix="$2" @@ -119,12 +123,28 @@ comment_out_line "$atten_input_file" "$atten_output_file" if [ "$SPECT" -eq 0 ]; then if [ "$TOF" -eq 0 ]; then - : ${view_mash:=1} - : ${span:=2} - : ${max_rd:=3} - echo "=== create template sinogram (DSTE with view_mash=${view_mash}, max_ring_diff=${max_rd})" - template_sino=my_DSTE_3D_vm${view_mash}_span${span}_rd${max_rd}_template.hs - cat > my_input.txt < my_input.txt < my_input.txt <1?argv[1]:""); + + + return reconstruction_object.reconstruct() == Succeeded::yes ? + EXIT_SUCCESS : EXIT_FAILURE; +} + diff --git a/src/analytic/DFM3D/DFM3DReconstruction.cxx b/src/analytic/DFM3D/DFM3DReconstruction.cxx new file mode 100644 index 000000000..978963d2a --- /dev/null +++ b/src/analytic/DFM3D/DFM3DReconstruction.cxx @@ -0,0 +1,629 @@ +/* + Copyright (C) 2024 University College London + + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup analytic + \brief Implementation of class stir::DFM3DReconstruction + + \author Dimitra Kyriakopoulou + \author Kris Thielemans +*/ + +#include "stir/analytic/DFM3D/DFM3DReconstruction.h" +#include "stir/VoxelsOnCartesianGrid.h" +#include "stir/RelatedViewgrams.h" +#include "stir/recon_buildblock/BackProjectorByBinUsingInterpolation.h" +#include "stir/ProjDataInfoCylindricalArcCorr.h" +#include "stir/ArcCorrection.h" +#include "stir/SSRB.h" +#include "stir/ProjDataInMemory.h" +#include "stir/Bin.h" +#include "stir/round.h" +#include "stir/display.h" +#include "stir/recon_buildblock/ArtificialScanner3D.h" +#include "stir/recon_buildblock/missing_data/MissingDataReprojection3D.h" +#include +#include "stir/IO/interfile.h" + +#include "stir/Sinogram.h" +#include "stir/Viewgram.h" +#include +#include +#include "stir/numerics/fourier.h" +#include +#include "stir/info.h" + +#ifdef STIR_OPENMP +#include +#endif +#include "stir/num_threads.h" + +#include "stir/Scanner.h" + +#include +#include +#include // For M_PI and other math functions +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +START_NAMESPACE_STIR + +const char * const +DFM3DReconstruction::registered_name = + "DFM3D"; + + +void +DFM3DReconstruction:: +set_defaults() +{ + base_type::set_defaults(); + display_level=0; // no display + num_segments_to_combine = -1; + noise_filter = 0; +} + +void +DFM3DReconstruction::initialise_keymap() +{ + base_type::initialise_keymap(); + + parser.add_start_key("DFM3DParameters"); + parser.add_stop_key("End"); + parser.add_key("num_segments_to_combine with SSRB", &num_segments_to_combine); + parser.add_key("Display level",&display_level); + parser.add_key("noise filter",&noise_filter); +} + +void +DFM3DReconstruction:: +ask_parameters() +{ + + base_type::ask_parameters(); + + num_segments_to_combine = ask_num("num_segments_to_combine (must be odd)",-1,101,-1); + display_level = ask_num("Which images would you like to display \n\t(0: None, 1: Final, 2: filtered viewgrams) ? ", 0,2,0); + noise_filter = ask_num(" Noise filter (0 to disable, 1 to enable)",0.,1., 0.); + } + +bool DFM3DReconstruction::post_processing() +{ + return base_type::post_processing(); +} + +Succeeded +DFM3DReconstruction:: +set_up(shared_ptr const& target_data_sptr) +{ + if (base_type::set_up(target_data_sptr) == Succeeded::no) + return Succeeded::no; + if (noise_filter != 1 && noise_filter != 0) + warning("Noise filter has to be either 0 or 1 but is %g\n", noise_filter); + + if (num_segments_to_combine>=0 && num_segments_to_combine%2==0) + error("num_segments_to_combine has to be odd (or -1), but is %d", num_segments_to_combine); + + if (num_segments_to_combine==-1) + { + const shared_ptr proj_data_info_cyl_sptr = + dynamic_pointer_cast(proj_data_ptr->get_proj_data_info_sptr()); + + if (is_null_ptr(proj_data_info_cyl_sptr)) + num_segments_to_combine = 1; //cannot SSRB non-cylindrical data yet + else + { + if (proj_data_info_cyl_sptr->get_min_ring_difference(0) != + proj_data_info_cyl_sptr->get_max_ring_difference(0) + || + proj_data_info_cyl_sptr->get_num_segments()==1) + num_segments_to_combine = 1; + else + num_segments_to_combine = 3; + } + } + + + return Succeeded::yes; +} + +std::string DFM3DReconstruction::method_info() const +{ + return "DFM3D"; +} + +DFM3DReconstruction:: +DFM3DReconstruction(const std::string& parameter_filename) +{ + initialise(parameter_filename); + //std::cerr<& proj_data_ptr_v, + const int noise_filter_v, + const int num_segments_to_combine_v +) +{ + set_defaults(); + + noise_filter = noise_filter_v; + num_segments_to_combine = num_segments_to_combine_v; + proj_data_ptr = proj_data_ptr_v; + // have to check here because we're not parsing + //if (post_processing_only_DFM3D_parameters() == true) + // error("DFM3D: Wrong parameter values. Aborting\n"); +} + +Succeeded +DFM3DReconstruction:: +actual_reconstruct(shared_ptr > const & density_ptr) +{ + + // perform SSRB +/* if (num_segments_to_combine>1) + { + const ProjDataInfoCylindrical& proj_data_info_cyl = + dynamic_cast + (*proj_data_ptr->get_proj_data_info_sptr()); + + // full_log << "SSRB combining " << num_segments_to_combine + // << " segments in input file to a new segment 0\n" << std::endl; + + shared_ptr + ssrb_info_sptr(SSRB(proj_data_info_cyl, + num_segments_to_combine, + 1, 0, + (num_segments_to_combine-1)/2 )); + shared_ptr + proj_data_to_GRD_ptr(new ProjDataInMemory (proj_data_ptr->get_exam_info_sptr(), ssrb_info_sptr)); + SSRB(*proj_data_to_GRD_ptr, *proj_data_ptr); + proj_data_ptr = proj_data_to_GRD_ptr; + } + else + { + // just use the proj_data_ptr we have already + }*/ + + // check if segment 0 has direct sinograms + { + const float tan_theta = proj_data_ptr->get_proj_data_info_sptr()->get_tantheta(Bin(0,0,0,0)); + if(fabs(tan_theta ) > 1.E-4) + { + warning("DFM3D: segment 0 has non-zero tan(theta) %g", tan_theta); + return Succeeded::no; + } + } + + // unused warning + float tangential_sampling; + // TODO make next type shared_ptr once we moved to boost::shared_ptr + // will enable us to get rid of a few of the ugly lines related to tangential_sampling below + shared_ptr arc_corrected_proj_data_info_sptr; + +/* // arc-correction if necessary + ArcCorrection arc_correction; + bool do_arc_correction = false; + if (!is_null_ptr(dynamic_pointer_cast + (proj_data_ptr->get_proj_data_info_sptr()))) + { + // it's already arc-corrected + arc_corrected_proj_data_info_sptr = + proj_data_ptr->get_proj_data_info_sptr()->create_shared_clone(); + tangential_sampling = + dynamic_cast + (*proj_data_ptr->get_proj_data_info_sptr()).get_tangential_sampling(); + } + else + { + // TODO arc-correct to voxel_size + if (arc_correction.set_up(proj_data_ptr->get_proj_data_info_sptr()->create_shared_clone()) == + Succeeded::no) + return Succeeded::no; + do_arc_correction = true; + // TODO full_log + warning("DFM3D will arc-correct data first"); + arc_corrected_proj_data_info_sptr = + arc_correction.get_arc_corrected_proj_data_info_sptr(); + tangential_sampling = + arc_correction.get_arc_corrected_proj_data_info().get_tangential_sampling(); + }*/ + //ProjDataInterfile ramp_filtered_proj_data(arc_corrected_proj_data_info_sptr,"ramp_filtered"); + + + VoxelsOnCartesianGrid& image = + dynamic_cast&>(*density_ptr); + + density_ptr->fill(0); + + const Scanner* scanner = proj_data_ptr->get_proj_data_info_sptr()->get_scanner_ptr(); + + Sinogram sino = proj_data_ptr->get_empty_sinogram(0,0); + const int spA = sino.get_num_tangential_poss(); + const int sth = sino.get_num_views(); + const ArtificialScanner3DLayout art_scanner_layout = + create_default_artificial_scanner3d_layout(*proj_data_ptr->get_proj_data_info_sptr()); + const int sphi = static_cast(art_scanner_layout.segment_numbers.size()); + if (sphi != proj_data_ptr->get_num_segments()) + error("DFM3D: artificial scanner layout mismatch (layout=%d data=%d)", sphi, proj_data_ptr->get_num_segments()); + const std::vector& segment_numbers = art_scanner_layout.segment_numbers; + const std::vector& axial = art_scanner_layout.target_axial_counts; + const int seg0 = art_scanner_layout.centre_index; + const int saA = art_scanner_layout.reference_axial_count; + const float dring = 2.0*scanner->get_ring_spacing()/(2.0*scanner->get_inner_ring_radius()); + const float dpA = 2.0/(spA-1); + const float dth = M_PI/sth; + + // pad to the next power of 2 for FFT + const int sp = 4*pow(2,ceil(log2(spA))); + const int sa = 2*pow(2,ceil(log2(saA))); + std::cout << "Tangential poss is " << spA << ", FFT length is " << sp << std::endl; + std::cout << "Axial poss is " << saA << ", FFT length is " << sa << std::endl; + + const float dp = 2.0/(sp-1); + + // IDWa parameters + const int L = 2; + const float p = 1.0; + // gridding + std::vector pn(sp,0.0f), an(sa,0.0f), thn(sth,0.0f), phin(sphi,0.0f), xn(sp,0.0f), yn(sp,0.0f), zn(sa,0.0f); + std::vector pnA(spA,0.0f), anA(saA,0.0f); + // angles + for(int ith=0; ith > vg(span); + IndexRange3D span1(sa,sp,sp); + Array< 3, std::complex > fg1(span1); + IndexRange3D span1A(saA,spA,spA); //p + Array< 3, std::complex > fg1A(span1A);//p + +std::cerr << "Version 240715" << std::endl; +if(image.get_x_size() != sino.get_num_tangential_poss()+((sino.get_num_tangential_poss()+1)%2)) { std::cerr << "Will interpolate" << std::endl; } + const int NP = sp*sphi*sth; + +std::vector> xt1(sa, std::vector(NP,0.0f)); +std::vector> yt1(sa, std::vector(NP,0.0f)); +std::vector>> ft1(sa, std::vector>(NP,.0f)); + // viewgr contains measured data plus artificial missing-data bins + MissingDataSinogram4D viewgr(sphi, sth, spA, axial); + embed_measured_viewgrams_into_missing_data_sinogram(viewgr, *proj_data_ptr, art_scanner_layout); + + // Calculate initial estimate + std::cout << std::endl << "Calculating initial image estimation... " << std::endl; + // -- calculate 2d fft of projected data + { + int iphi = seg0, i = 0; + + for(int ith = 0; ith hanningA(sa,0.0f), hanningP1(sp,0.0f), hanningP2(sp,0.0f); + for (int i = 0; i < sa; ++i) + hanningA[i] = 0.5 * (1 - cos(2 * M_PI * i / (sa - 1))); + for (int j = 0; j < sp; ++j) + hanningP1[j] = 0.5 * (1 - cos(2 * M_PI * j / (sp - 1))); + for (int k = 0; k < sp; ++k) + hanningP2[k] = 0.5 * (1 - cos(2 * M_PI * k / (sp - 1))); + // Apply the Hanning filter to the data in all three dimensions + for (int i = 0; i < sa; ++i) { + for (int j = 0; j < sp; ++j) { + for (int k = 0; k < sp; ++k) { + // Multiplicative application of the Hanning filter + fg1[i][j][k] *= hanningA[i] * hanningP1[j] * hanningP2[k]; + } + } + } + } + + + + // -- calculate inverse 3d fft + fftshift(fg1,sa,sp,sp); + inverse_fourier(fg1); + fftshift(fg1,sa,sp,sp); + + + // Reproject //p + fg1A.fill(.0f); + for(int iz=0; izget_viewgram(ith,seg); + + vg.fill(.0f); + // copy data + for(int ip=0; ip xn1(sx,0.0f); + std::vector yn1(sy,0.0f); + float dx1 = 2./(sx-1), dy1 = 2./(sy-1); + float dx = 2./(spA-1), dy = 2./(spA-1); + float val; + + for(int ix=0; ix0) + display(image, image.find_max(), "DFM3D image"); + + return Succeeded::yes; +} + + + +template +void +DFM3DReconstruction::fftshift(Array< 1 , T >& a, int size) +{ + T temp=0; + for(int i=0; i +void +DFM3DReconstruction::fftshift(Array< 2 , std::complex< T > >& a, int M, int N) +{ + std::complex temp; + for(int i=0; i +void +DFM3DReconstruction::fftshift(Array< 3 , std::complex< T > >& a, int M, int N, int K) +{ + std::complex temp; + for(int i=0; i > DFM3DReconstruction::IDWa(const std::vector& xt1, const std::vector& yt1, const std::vector>& ft1, int N, const std::vector& xn, const std::vector& yn, int L, float p, int sp) +{ + IndexRange2D span(sp,sp); + Array< 2, std::complex > fg1(span); + +std::vector> fg1h(sp, std::vector(sp, 0.0f)); +for(int i=0; i +#include "stir/shared_ptr.h" +#include "stir/Array_complex_numbers.h" + + +START_NAMESPACE_STIR + +template class DiscretisedDensity; +class Succeeded; +class ProjData; + +class DFM3DReconstruction : + public + RegisteredParsingObject< + DFM3DReconstruction, + Reconstruction < DiscretisedDensity < 3,float> >, + AnalyticReconstruction + > +{ + //typedef AnalyticReconstruction base_type; + typedef + RegisteredParsingObject< + DFM3DReconstruction, + Reconstruction < DiscretisedDensity < 3,float> >, + AnalyticReconstruction + > base_type; + typedef DiscretisedDensity < 3,float> TargetT; +public: + //! Name which will be used when parsing a ProjectorByBinPair object + static const char * const registered_name; + + //! Default constructor (calls set_defaults()) + DFM3DReconstruction (); + /*! + \brief Constructor, initialises everything from parameter file, or (when + parameter_filename == "") by calling ask_parameters(). + */ + explicit + DFM3DReconstruction(const std::string& parameter_filename); + + DFM3DReconstruction(const shared_ptr&, + const int noise_filter=0, + const int num_segments_to_combine=-1 + ); + + virtual std::string method_info() const; + + virtual void ask_parameters(); + + virtual Succeeded set_up(shared_ptr const& target_data_sptr); + + protected: // make parameters protected such that doc shows always up in doxygen + // parameters used for parsing + + // noise filter + int noise_filter; + //! number of segments to combine (with SSRB) before starting 2D reconstruction + /*! if -1, a value is chosen depending on the axial compression. + If there is no axial compression, num_segments_to_combine is + effectively set to 3, otherwise it is set to 1. + \see SSRB + */ + int num_segments_to_combine; + //! potentially display data + /*! allowed values: \c display_level=0 (no display), 1 (only final image), + 2 (filtered-viewgrams). Defaults to 0. + */ + int display_level; + private: + Succeeded actual_reconstruct(shared_ptr > const & target_image_ptr); + + + virtual void set_defaults(); + virtual void initialise_keymap(); + virtual bool post_processing(); + bool post_processing_only_DFM3D_parameters(); + +void process_chunk(int start, int end, int sp, int sphi, int sth, int sa); + template void fftshift(Array< 1 , T >& a, int size); + template void fftshift(Array< 2 , std::complex< T > >& a, int M, int N); + template void fftshift(Array< 3 , std::complex< T > >& a, int M, int N, int K); +Array< 2, std::complex > IDWa(const std::vector& xt1, const std::vector& yt1, const std::vector>& ft1, int N, const std::vector& xn, const std::vector& yn, int L, float p, int sp); +}; + + + + +END_NAMESPACE_STIR + + +#endif diff --git a/src/include/stir/recon_buildblock/ArtificialScanner3D.h b/src/include/stir/recon_buildblock/ArtificialScanner3D.h new file mode 100644 index 000000000..2f9bd1b37 --- /dev/null +++ b/src/include/stir/recon_buildblock/ArtificialScanner3D.h @@ -0,0 +1,56 @@ +/* + Copyright (C) 2026 University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup recon_buildblock + + \brief Utilities to derive a reusable artificial scanner axial layout for 3D sinograms. + +*/ + +#ifndef __stir_recon_buildblock_ArtificialScanner3D_H__ +#define __stir_recon_buildblock_ArtificialScanner3D_H__ + +#include "stir/common.h" +#include + +START_NAMESPACE_STIR + +class ProjDataInfo; + +/*! + \brief Compact per-segment layout used by artificial-scanner based 3D algorithms. + + The vectors are aligned by segment index \c iphi. +*/ +struct ArtificialScanner3DLayout +{ + std::vector segment_numbers; + std::vector measured_axial_counts; + std::vector target_axial_counts; + std::vector measured_offsets; + int centre_index = 0; + int centre_segment_num = 0; + int reference_axial_count = 0; +}; + +/*! + \brief Construct a default artificial scanner layout from projection-data geometry. + + The target axial count for each segment is: + \code + max(measured_axial_count, reference_axial_count + |seg - centre_segment|) + \endcode +*/ +ArtificialScanner3DLayout +create_default_artificial_scanner3d_layout(const ProjDataInfo& proj_data_info); + +END_NAMESPACE_STIR + +#endif diff --git a/src/include/stir/recon_buildblock/missing_data/MissingDataReprojection3D.h b/src/include/stir/recon_buildblock/missing_data/MissingDataReprojection3D.h new file mode 100644 index 000000000..176979cfc --- /dev/null +++ b/src/include/stir/recon_buildblock/missing_data/MissingDataReprojection3D.h @@ -0,0 +1,90 @@ +/* + Copyright (C) 2026 University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup recon_buildblock + + \brief Standalone utilities for DFM-style 3D missing-data filling. + +*/ + +#ifndef __stir_recon_buildblock_missing_data_MissingDataReprojection3D_H__ +#define __stir_recon_buildblock_missing_data_MissingDataReprojection3D_H__ + +#include "stir/Array_complex_numbers.h" +#include "stir/common.h" +#include +#include + +START_NAMESPACE_STIR + +class ProjData; +struct ArtificialScanner3DLayout; + +/*! + \brief Dense 4D storage for per-segment sinograms, including artificial missing-data bins. + + Layout indexes are \c (segment_index, view, axial, tangential). +*/ +class MissingDataSinogram4D +{ +public: + MissingDataSinogram4D() = default; + MissingDataSinogram4D(int num_segments, int num_views, int num_tangential_poss, const std::vector& axial_counts); + + void + resize(int num_segments, int num_views, int num_tangential_poss, const std::vector& axial_counts); + + inline int get_num_segments() const { return static_cast(this->axial_counts_.size()); } + inline int get_num_views() const { return this->num_views_; } + inline int get_num_tangential_poss() const { return this->num_tangential_poss_; } + inline int get_num_axial_poss(const int segment_index) const { return this->axial_counts_[segment_index]; } + + float& + at(int segment_index, int view, int axial, int tangential); + const float& + at(int segment_index, int view, int axial, int tangential) const; + +private: + std::size_t + compute_offset(int segment_index, int view, int axial, int tangential) const; + + std::vector axial_counts_; + std::vector segment_offsets_; + int num_views_ = 0; + int num_tangential_poss_ = 0; + std::vector data_; +}; + +/*! + \brief Copy measured projection data into a missing-data container according to artificial-scanner layout. +*/ +void +embed_measured_viewgrams_into_missing_data_sinogram( + MissingDataSinogram4D& destination, const ProjData& proj_data, const ArtificialScanner3DLayout& layout); + +/*! + \brief Fill only missing bins with DFM-style trilinear reprojection from \c image. + + Measured bins are left untouched, defined by \c layout.measured_offsets and + \c layout.measured_axial_counts. +*/ +void +fill_missing_data_by_trilinear_reprojection(MissingDataSinogram4D& destination, + const Array<3, std::complex>& image, + const std::vector& pn, + const std::vector& an, + const std::vector& thn, + const std::vector& phin, + const ArtificialScanner3DLayout& layout, + std::ostream* log_stream = nullptr); + +END_NAMESPACE_STIR + +#endif diff --git a/src/recon_buildblock/ArtificialScanner3D.cxx b/src/recon_buildblock/ArtificialScanner3D.cxx new file mode 100644 index 000000000..b3b8fb187 --- /dev/null +++ b/src/recon_buildblock/ArtificialScanner3D.cxx @@ -0,0 +1,83 @@ +/* + Copyright (C) 2026 University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup recon_buildblock + + \brief Implementation of artificial scanner layout utilities for 3D sinograms. + +*/ + +#include "stir/recon_buildblock/ArtificialScanner3D.h" +#include "stir/ProjDataInfo.h" +#include "stir/error.h" +#include "stir/warning.h" +#include +#include + +START_NAMESPACE_STIR + +ArtificialScanner3DLayout +create_default_artificial_scanner3d_layout(const ProjDataInfo& proj_data_info) +{ + ArtificialScanner3DLayout layout; + + const int min_segment_num = proj_data_info.get_min_segment_num(); + const int max_segment_num = proj_data_info.get_max_segment_num(); + if (max_segment_num < min_segment_num) + error("ArtificialScanner3D: invalid segment range (min=%d max=%d)", min_segment_num, max_segment_num); + + const int num_segments = max_segment_num - min_segment_num + 1; + layout.segment_numbers.resize(num_segments); + for (int i = 0; i < num_segments; ++i) + layout.segment_numbers[i] = min_segment_num + i; + + layout.centre_index = num_segments / 2; + if (min_segment_num <= 0 && 0 <= max_segment_num) + layout.centre_index = -min_segment_num; + if (layout.centre_index < 0 || layout.centre_index >= num_segments) + error("ArtificialScanner3D: centre index out of range (%d for %d segments)", layout.centre_index, num_segments); + layout.centre_segment_num = layout.segment_numbers[layout.centre_index]; + + if (min_segment_num <= 0 && 0 <= max_segment_num) + layout.reference_axial_count = proj_data_info.get_num_axial_poss(0); + if (layout.reference_axial_count <= 0) + { + for (int i = 0; i < num_segments; ++i) + layout.reference_axial_count + = std::max(layout.reference_axial_count, proj_data_info.get_num_axial_poss(layout.segment_numbers[i])); + warning("ArtificialScanner3D: segment 0 unavailable; using max axial count %d as reference", + layout.reference_axial_count); + } + if (layout.reference_axial_count <= 0) + error("ArtificialScanner3D: failed to derive a positive reference axial count"); + + layout.measured_axial_counts.resize(num_segments); + layout.target_axial_counts.resize(num_segments); + layout.measured_offsets.resize(num_segments); + + for (int i = 0; i < num_segments; ++i) + { + const int segment_num = layout.segment_numbers[i]; + const int measured_count = proj_data_info.get_num_axial_poss(segment_num); + if (measured_count <= 0) + error("ArtificialScanner3D: non-positive axial count %d at segment %d", measured_count, segment_num); + + const int ring_diff_from_centre = std::abs(segment_num - layout.centre_segment_num); + const int target_count = std::max(measured_count, layout.reference_axial_count + ring_diff_from_centre); + + layout.measured_axial_counts[i] = measured_count; + layout.target_axial_counts[i] = target_count; + layout.measured_offsets[i] = (target_count - measured_count) / 2; + } + + return layout; +} + +END_NAMESPACE_STIR diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index b0d6466df..a44889aeb 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -14,6 +14,8 @@ set(dir recon_buildblock) set (dir_LIB_SOURCES ${dir}_LIB_SOURCES) set(${dir_LIB_SOURCES} + ArtificialScanner3D.cxx + missing_data/MissingDataReprojection3D.cxx ForwardProjectorByBin.cxx BackProjectorByBin.cxx ProjectorByBinPair.cxx diff --git a/src/recon_buildblock/missing_data/MissingDataReprojection3D.cxx b/src/recon_buildblock/missing_data/MissingDataReprojection3D.cxx new file mode 100644 index 000000000..e71269ef3 --- /dev/null +++ b/src/recon_buildblock/missing_data/MissingDataReprojection3D.cxx @@ -0,0 +1,250 @@ +/* + Copyright (C) 2026 University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup recon_buildblock + + \brief Implementation of standalone DFM-style 3D missing-data filling utilities. + +*/ + +#include "stir/recon_buildblock/missing_data/MissingDataReprojection3D.h" +#include "stir/recon_buildblock/ArtificialScanner3D.h" +#include "stir/ProjData.h" +#include "stir/ProjDataInfo.h" +#include "stir/Viewgram.h" +#include "stir/error.h" +#include +#include +#include + +START_NAMESPACE_STIR + +MissingDataSinogram4D::MissingDataSinogram4D( + int num_segments, int num_views, int num_tangential_poss, const std::vector& axial_counts) +{ + this->resize(num_segments, num_views, num_tangential_poss, axial_counts); +} + +void +MissingDataSinogram4D::resize( + const int num_segments, const int num_views, const int num_tangential_poss, const std::vector& axial_counts) +{ + if (num_segments <= 0) + error("MissingDataSinogram4D: num_segments must be positive (%d)", num_segments); + if (num_views <= 0) + error("MissingDataSinogram4D: num_views must be positive (%d)", num_views); + if (num_tangential_poss <= 0) + error("MissingDataSinogram4D: num_tangential_poss must be positive (%d)", num_tangential_poss); + if (static_cast(axial_counts.size()) != num_segments) + error("MissingDataSinogram4D: axial_counts size mismatch (%d vs %d)", static_cast(axial_counts.size()), num_segments); + + this->axial_counts_ = axial_counts; + this->num_views_ = num_views; + this->num_tangential_poss_ = num_tangential_poss; + + this->segment_offsets_.assign(static_cast(num_segments) + 1U, 0U); + for (int iphi = 0; iphi < num_segments; ++iphi) + { + if (this->axial_counts_[iphi] <= 0) + error("MissingDataSinogram4D: non-positive axial count %d at segment index %d", this->axial_counts_[iphi], iphi); + + const std::size_t segment_size = static_cast(num_views) * static_cast(this->axial_counts_[iphi]) + * static_cast(num_tangential_poss); + this->segment_offsets_[static_cast(iphi) + 1U] + = this->segment_offsets_[static_cast(iphi)] + segment_size; + } + + this->data_.assign(this->segment_offsets_.back(), 0.F); +} + +std::size_t +MissingDataSinogram4D::compute_offset(const int segment_index, const int view, const int axial, const int tangential) const +{ + if (segment_index < 0 || segment_index >= this->get_num_segments()) + error("MissingDataSinogram4D: segment index out of range (%d)", segment_index); + if (view < 0 || view >= this->num_views_) + error("MissingDataSinogram4D: view index out of range (%d)", view); + if (axial < 0 || axial >= this->axial_counts_[segment_index]) + error("MissingDataSinogram4D: axial index out of range (%d)", axial); + if (tangential < 0 || tangential >= this->num_tangential_poss_) + error("MissingDataSinogram4D: tangential index out of range (%d)", tangential); + + return this->segment_offsets_[static_cast(segment_index)] + + (static_cast(view) * static_cast(this->axial_counts_[segment_index]) + + static_cast(axial)) + * static_cast(this->num_tangential_poss_) + + static_cast(tangential); +} + +float& +MissingDataSinogram4D::at(const int segment_index, const int view, const int axial, const int tangential) +{ + return this->data_[this->compute_offset(segment_index, view, axial, tangential)]; +} + +const float& +MissingDataSinogram4D::at(const int segment_index, const int view, const int axial, const int tangential) const +{ + return this->data_[this->compute_offset(segment_index, view, axial, tangential)]; +} + +void +embed_measured_viewgrams_into_missing_data_sinogram( + MissingDataSinogram4D& destination, const ProjData& proj_data, const ArtificialScanner3DLayout& layout) +{ + const int sphi = static_cast(layout.segment_numbers.size()); + if (destination.get_num_segments() != sphi) + error("embed_measured_viewgrams_into_missing_data_sinogram: segment mismatch (%d vs %d)", + destination.get_num_segments(), + sphi); + if (static_cast(layout.measured_axial_counts.size()) != sphi || static_cast(layout.measured_offsets.size()) != sphi) + error("embed_measured_viewgrams_into_missing_data_sinogram: inconsistent layout vectors"); + + const int min_tangential = proj_data.get_proj_data_info_sptr()->get_min_tangential_pos_num(); + + for (int iphi = 0; iphi < sphi; ++iphi) + { + const int seg = layout.segment_numbers[iphi]; + const int measured_axial = layout.measured_axial_counts[iphi]; + const int axial_offset = layout.measured_offsets[iphi]; + + for (int ith = 0; ith < destination.get_num_views(); ++ith) + { + const Viewgram view = proj_data.get_viewgram(ith, seg); + const int min_axial = view.get_min_axial_pos_num(); + + for (int ia = 0; ia < measured_axial; ++ia) + { + for (int ip = 0; ip < destination.get_num_tangential_poss(); ++ip) + { + destination.at(iphi, ith, ia + axial_offset, ip) = view[ia + min_axial][ip + min_tangential]; + } + } + } + } +} + +void +fill_missing_data_by_trilinear_reprojection(MissingDataSinogram4D& destination, + const Array<3, std::complex>& image, + const std::vector& pn, + const std::vector& an, + const std::vector& thn, + const std::vector& phin, + const ArtificialScanner3DLayout& layout, + std::ostream* log_stream) +{ + const int sphi = static_cast(layout.segment_numbers.size()); + if (destination.get_num_segments() != sphi) + error("fill_missing_data_by_trilinear_reprojection: segment mismatch (%d vs %d)", destination.get_num_segments(), sphi); + if (static_cast(layout.target_axial_counts.size()) != sphi || static_cast(layout.measured_offsets.size()) != sphi + || static_cast(layout.measured_axial_counts.size()) != sphi) + error("fill_missing_data_by_trilinear_reprojection: inconsistent layout vectors"); + + const int sp = static_cast(pn.size()); + const int sa = static_cast(an.size()); + const int sth = static_cast(thn.size()); + if (sp < 2 || sa < 2) + error("fill_missing_data_by_trilinear_reprojection: pn/an sizes must be >=2 (sp=%d sa=%d)", sp, sa); + if (destination.get_num_tangential_poss() != sp) + error("fill_missing_data_by_trilinear_reprojection: tangential size mismatch (%d vs %d)", + destination.get_num_tangential_poss(), + sp); + if (destination.get_num_views() != sth) + error("fill_missing_data_by_trilinear_reprojection: view size mismatch (%d vs %d)", destination.get_num_views(), sth); + if (static_cast(phin.size()) != sphi) + error("fill_missing_data_by_trilinear_reprojection: phin size mismatch (%d vs %d)", static_cast(phin.size()), sphi); + + const float dp = pn[1] - pn[0]; + const float da = an[1] - an[0]; + if (dp == 0.F || da == 0.F) + error("fill_missing_data_by_trilinear_reprojection: zero sampling step (dp=%g da=%g)", dp, da); + + if (log_stream) + *log_stream << "Reprojecting... [segment_num(num_axial_poss)]" << std::endl; + + for (int iphi = 0; iphi < sphi; ++iphi) + { + const int seg = layout.segment_numbers[iphi]; + const int sa1 = layout.target_axial_counts[iphi]; + const int known_begin = layout.measured_offsets[iphi]; + const int known_end = known_begin + layout.measured_axial_counts[iphi]; + + if (destination.get_num_axial_poss(iphi) != sa1) + error("fill_missing_data_by_trilinear_reprojection: axial size mismatch at seg index %d (%d vs %d)", + iphi, + destination.get_num_axial_poss(iphi), + sa1); + + if (log_stream) + *log_stream << seg << "(" << layout.measured_axial_counts[iphi] << "->" << sa1 << "), " << std::flush; + + std::vector an1(sa1, 0.F); + for (int ia = 0; ia < sa1; ++ia) + an1[ia] = std::cos(phin[iphi]) * (-(sa1 - 1.0F) / 2.0F * da + ia * da); + + for (int ith = 0; ith < sth; ++ith) + { + for (int ia = 0; ia < sa1; ++ia) + { + if (ia >= known_begin && ia < known_end) + continue; + + for (int ip = 0; ip < sp; ++ip) + { + for (int is = 0; is < sp; ++is) + { + const float x = pn[is] * std::cos(thn[ith]) * std::cos(phin[iphi]) - pn[ip] * std::sin(thn[ith]) + - an1[ia] * std::cos(thn[ith]) * std::sin(phin[iphi]); + const float y = pn[is] * std::sin(thn[ith]) * std::cos(phin[iphi]) + pn[ip] * std::cos(thn[ith]) + - an1[ia] * std::sin(thn[ith]) * std::sin(phin[iphi]); + const float z = pn[is] * std::sin(phin[iphi]) + an1[ia] * std::cos(phin[iphi]); + + if (x < -1.F || x > 1.F || y < -1.F || y > 1.F || z < an.front() || z > an.back()) + continue; + + int i = static_cast(std::floor((sp - 1) * (x + 1.0F) / 2.0F)); + if (x > 1.0F - 2.0e-6F) + i = sp - 2; + + int j = static_cast(std::floor((sp - 1) * (y + 1.0F) / 2.0F)); + if (y > 1.0F - 2.0e-6F) + j = sp - 2; + + int k = static_cast(std::floor((sa - 1) * (z + an.back()) / (an.back() - an.front()))); + if (z > an.back() - 2.0e-6F) + k = sa - 2; + + const float xd = (x - pn[i]) / dp; + const float yd = (y - pn[j]) / dp; + const float zd = (z - an[k]) / da; + + const float f00 = std::real(image[k][i][j]) * (1.0F - xd) + std::real(image[k][i + 1][j]) * xd; + const float f01 = std::real(image[k + 1][i][j]) * (1.0F - xd) + std::real(image[k + 1][i + 1][j]) * xd; + const float f10 = std::real(image[k][i][j + 1]) * (1.0F - xd) + std::real(image[k][i + 1][j + 1]) * xd; + const float f11 + = std::real(image[k + 1][i][j + 1]) * (1.0F - xd) + std::real(image[k + 1][i + 1][j + 1]) * xd; + + const float f0 = f00 * (1.0F - yd) + f10 * yd; + const float f1 = f01 * (1.0F - yd) + f11 * yd; + const float f = f0 * (1.0F - zd) + f1 * zd; + + float& out_bin = destination.at(iphi, ith, ia, ip); + out_bin += f; + if (is == 0 || is == sp - 1) + out_bin -= 0.5F * f; + } + } + } + } + } +} + +END_NAMESPACE_STIR diff --git a/src/recon_buildblock/recon_buildblock_registries.cxx b/src/recon_buildblock/recon_buildblock_registries.cxx index c6c35a6b4..a64af31cf 100644 --- a/src/recon_buildblock/recon_buildblock_registries.cxx +++ b/src/recon_buildblock/recon_buildblock_registries.cxx @@ -57,6 +57,7 @@ #include "stir/analytic/FBP2D/FBP2DReconstruction.h" #include "stir/analytic/FBP3DRP/FBP3DRPReconstruction.h" +#include "stir/analytic/DFM3D/DFM3DReconstruction.h" #include "stir/OSMAPOSL/OSMAPOSLReconstruction.h" #include "stir/KOSMAPOSL/KOSMAPOSLReconstruction.h" @@ -133,6 +134,7 @@ static PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion>::RegisterIt dummy603; static KOSMAPOSLReconstruction>::RegisterIt dummyK;