From b06adca995dcacfd6318184cd42098d89dd0e75b Mon Sep 17 00:00:00 2001 From: Daniel Deidda Date: Wed, 4 Mar 2026 17:01:02 +0000 Subject: [PATCH 1/9] SPECTGPU:starting point; created SPECTGPU_projector folder and scheleton classe, atm a lot of copy from niftyGPU --- .../BackProjectorByBinSPECTGPU.h | 96 +++ .../ForwardProjectorByBinSPECTGPU.h | 108 ++++ .../ProjectorByBinPairUsingSPECTGPU.h | 64 ++ .../SPECTGPU_projector/SPECTGPUHelper.h | 148 +++++ src/recon_buildblock/CMakeLists.txt | 4 + .../BackProjectorByBinSPECTGPU.cxx | 130 ++++ .../ForwardProjectorByBinSPECTGPU.cxx | 132 ++++ .../ProjectorByBinPairUsingSPECTGPU.cxx | 98 +++ .../SPECTGPU_projector/SPECTGPUHelper.cxx | 609 ++++++++++++++++++ .../SPECTGPU_projector/SPECTGPUProjection.cu | 52 ++ .../SPECTGPURotateAndGaussianInterpolate.cu | 51 ++ 11 files changed, 1492 insertions(+) create mode 100644 src/include/stir/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.h create mode 100644 src/include/stir/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.h create mode 100644 src/include/stir/recon_buildblock/SPECTGPU_projector/ProjectorByBinPairUsingSPECTGPU.h create mode 100644 src/include/stir/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.h create mode 100644 src/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.cxx create mode 100644 src/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.cxx create mode 100644 src/recon_buildblock/SPECTGPU_projector/ProjectorByBinPairUsingSPECTGPU.cxx create mode 100644 src/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.cxx create mode 100644 src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu create mode 100644 src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu diff --git a/src/include/stir/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.h b/src/include/stir/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.h new file mode 100644 index 0000000000..023cc75cab --- /dev/null +++ b/src/include/stir/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.h @@ -0,0 +1,96 @@ +// +// +/*! + \file + \ingroup projection + \ingroup SPECTGPU + + \brief Back projection class using SPECTGPU's GPU implementation. + + \author Daniel Deidda + + \todo SPECTGPU limitations - + + \todo STIR wrapper limitations - +*/ +/* + Copyright (C) 2026, National Physical Laboratory + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +#ifndef __stir_gpu_BackProjectorByBinSPECTGPU_h__ +#define __stir_gpu_BackProjectorByBinSPECTGPU_h__ + +#include "stir/RegisteredParsingObject.h" +#include "stir/recon_buildblock/BackProjectorByBin.h" +#include "stir/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.h" + +START_NAMESPACE_STIR + +class DataSymmetriesForViewSegmentNumbers; + +/*! + \ingroup projection + \brief Class for SPECTGPU's GPU back projector +*/ +class BackProjectorByBinSPECTGPU : public RegisteredParsingObject +{ +public: + //! Name which will be used when parsing a BackProjectorByBin object + static const char* const registered_name; + + //! Default constructor calls reset_timers() + BackProjectorByBinSPECTGPU(); + + virtual ~BackProjectorByBinSPECTGPU(); + + /// Keymap + virtual void initialise_keymap(); + + //! Stores all necessary geometric info + /*! + If necessary, set_up() can be called more than once. + */ + virtual void set_up(const shared_ptr& proj_data_info_ptr, + const shared_ptr>& density_info_sptr // TODO should be Info only + ); + + /// Back project + void back_project(const ProjData&, int subset_num = 0, int num_subsets = 1); + + /// Get output + virtual void get_output(DiscretisedDensity<3, float>&) const; + + /*! \brief tell the back projector to start accumulating into a new target. + This function has to be called before any back-projection is initiated.*/ + virtual void start_accumulating_in_new_target(); + + /// Set verbosity + void set_verbosity(const bool verbosity) { _cuda_verbosity = verbosity; } + + /// Set use truncation - truncate before forward + /// projection and after back projection + void set_use_truncation(const bool use_truncation) { _use_truncation = use_truncation; } + +protected: + virtual void actual_back_project(const RelatedViewgrams&, + const int min_axial_pos_num, + const int max_axial_pos_num, + const int min_tangential_pos_num, + const int max_tangential_pos_num); + +private: + shared_ptr _symmetries_sptr; + SPECTGPUHelper _helper; + int _cuda_device; + bool _cuda_verbosity; + std::vector _np_sino; + bool _use_truncation; +}; + +END_NAMESPACE_STIR + +#endif // __stir_gpu_BackProjectorByBinSPECTGPU_h__ diff --git a/src/include/stir/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.h b/src/include/stir/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.h new file mode 100644 index 0000000000..93868fe953 --- /dev/null +++ b/src/include/stir/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.h @@ -0,0 +1,108 @@ +// +// + +#ifndef __stir_gpu_ForwardProjectorByBinSPECTGPU_h__ +#define __stir_gpu_ForwardProjectorByBinSPECTGPU_h__ +/*! + \file + \ingroup projection + \ingroup SPECTGPU + + \brief Forward projection class using SPECTGPU's GPU implementation. + + \author Daniel Deidda + + \todo SPECTGPU limitations - + + \todo STIR wrapper limitations - + +*/ +/* + Copyright (C) 2026, National Physical Laboratory + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#include "stir/RegisteredParsingObject.h" +#include "stir/recon_buildblock/ForwardProjectorByBin.h" +#include "stir/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.h" + +START_NAMESPACE_STIR + +class ProjDataInMemory; +class DataSymmetriesForViewSegmentNumbers; + +/*! + \ingroup projection + \brief Class for SPECTGPU's GPU forward projector. +*/ +class ForwardProjectorByBinSPECTGPU : public RegisteredParsingObject +{ +public: + //! Name which will be used when parsing a ForwardProjectorByBin object + static const char* const registered_name; + + //! Default constructor calls reset_timers() + // inline + ForwardProjectorByBinSPECTGPU(); + + /// Constructor + virtual ~ForwardProjectorByBinSPECTGPU(); + + /// Keymap + virtual void initialise_keymap(); + + //! Stores all necessary geometric info + /*! + If necessary, set_up() can be called more than once. + + Derived classes can assume that forward_project() will be called + with input corresponding to the arguments of the last call to set_up(). + + \warning there is currently no check on this. + \warning Derived classes have to call set_up from the base class. + */ + virtual void set_up(const shared_ptr& proj_data_info_ptr, + const shared_ptr>& density_info_sptr // TODO should be Info only + ); + + /// Set input + virtual void set_input(const DiscretisedDensity<3, float>&); + + /// Set verbosity + void set_verbosity(const bool verbosity) { _cuda_verbosity = verbosity; } + + /// Set use truncation - truncate before forward + /// projection and after back projection + void set_use_truncation(const bool use_truncation) { _use_truncation = use_truncation; } + +protected: + //! This virtual function has to be implemented by the derived class. + virtual void actual_forward_project(RelatedViewgrams&, + const DiscretisedDensity<3, float>&, + const int min_axial_pos_num, + const int max_axial_pos_num, + const int min_tangential_pos_num, + const int max_tangential_pos_num); + + virtual void actual_forward_project(RelatedViewgrams& viewgrams, + const int min_axial_pos_num, + const int max_axial_pos_num, + const int min_tangential_pos_num, + const int max_tangential_pos_num); + +private: + shared_ptr _symmetries_sptr; + shared_ptr _projected_data_sptr; + SPECTGPUHelper _helper; + int _cuda_device; + bool _cuda_verbosity; + bool _use_truncation; +}; + +END_NAMESPACE_STIR + +#endif // __stir_gpu_ForwardProjectorByBinSPECTGPU_h__ diff --git a/src/include/stir/recon_buildblock/SPECTGPU_projector/ProjectorByBinPairUsingSPECTGPU.h b/src/include/stir/recon_buildblock/SPECTGPU_projector/ProjectorByBinPairUsingSPECTGPU.h new file mode 100644 index 0000000000..2b0742dc62 --- /dev/null +++ b/src/include/stir/recon_buildblock/SPECTGPU_projector/ProjectorByBinPairUsingSPECTGPU.h @@ -0,0 +1,64 @@ +// +// +/* + Copyright (C) 2026, National Physical Laboratory + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup projection + \ingroup SPECTGPU + + \brief Declares class stir::ProjectorByBinPairUsingSPECTGPU + + \author Daniel Deidda + +*/ +#ifndef __stir_recon_buildblock_ProjectorByBinPairUsingSPECTGPU_h_ +#define __stir_recon_buildblock_ProjectorByBinPairUsingSPECTGPU_h_ + +#include "stir/RegisteredParsingObject.h" +#include "stir/recon_buildblock/ProjectorByBinPair.h" + +START_NAMESPACE_STIR + +class Succeeded; +/*! + \ingroup projection + \brief A projector pair based on SPECTGPU projectors +*/ +class ProjectorByBinPairUsingSPECTGPU + : public RegisteredParsingObject +{ +private: + typedef RegisteredParsingObject base_type; + +public: + //! Name which will be used when parsing a ProjectorByBinPair object + static const char* const registered_name; + + //! Default constructor + ProjectorByBinPairUsingSPECTGPU(); + + /// Set verbosity + void set_verbosity(const bool verbosity); + + /// Set use truncation - truncate before forward + /// projection and after back projection + void set_use_truncation(const bool use_truncation); + +private: + void set_defaults(); + void initialise_keymap(); + bool post_processing(); + bool _verbosity; + bool _use_truncation; +}; + +END_NAMESPACE_STIR + +#endif // __stir_recon_buildblock_ProjectorByBinPairUsingSPECTGPU_h_ diff --git a/src/include/stir/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.h b/src/include/stir/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.h new file mode 100644 index 0000000000..2d77df3ae1 --- /dev/null +++ b/src/include/stir/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.h @@ -0,0 +1,148 @@ +// +// + +#ifndef __stir_gpu_SPECTGPUHelper_h__ +#define __stir_gpu_SPECTGPUHelper_h__ +/*! + \file + \ingroup projection + \ingroup SPECTGPU + + \brief Helper class for SPECTGPU's GPU implementation. + + \author Daniel Deidda + + Helper class for SPECTGPU's GPU functionality. Wrapped + functionality includes projection. + + \todo SPECTGPU limitations - . + + \todo STIR wrapper limitations - +/* + Copyright (C) 2026, National Physical Laboratory + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#include "stir/shared_ptr.h" +#include "stir/Scanner.h" + +// Forward declarations +struct Cnst; +struct txLUTs; +struct axialLUT; + +START_NAMESPACE_STIR + +template +class DiscretisedDensity; +class ProjData; +template +class Viewgram; +template +class VoxelsOnCartesianGrid; + +/*! + \ingroup projection + \brief Helper class for the wrapped SPECTGPU projectors. +*/ +class SPECTGPUHelper +{ +public: + /// Default constructor + SPECTGPUHelper() + : _already_set_up(false), + _devid(0), + _att(-1), + _scanner_type(Scanner::Unknown_scanner) + {} + + /// Destructor + virtual ~SPECTGPUHelper(); + + /// Set CUDA device ID + void set_cuda_device_id(const int devid) { _devid = char(devid); } + + + /// Set emission (0) or transmission (1) - whether to exp{-result} for attenuation maps + void set_att(const char att) { _att = att; } + + /// Set verbosity level for CUDA output + void set_verbose(const bool verbose) { _verbose = verbose; } + + /// Set scanner type + void set_scanner_type(const Scanner::Type scanner_type) { _scanner_type = scanner_type; } + + /// Set up + void set_up(); + + /// Create SPECTGPU image + static std::vector create_SPECTGPU_image(); + + /// Create STIR image with template dim + static shared_ptr> create_stir_im(); + + /// Create SPECTGPU singram. Forward project into this + std::vector create_SPECTGPU_sinogram() const; + + /// Convert STIR image to SPECTGPU image + static void convert_image_stir_to_SPECTGPU(std::vector& np, const DiscretisedDensity<3, float>& stir); + + /// Convert SPECTGPU image to STIR image + static void convert_image_SPECTGPU_to_stir(DiscretisedDensity<3, float>& stir, const std::vector& np_vec); + + /// Convert STIR proj data to SPECTGPU proj data + void convert_proj_data_stir_to_SPECTGPU(std::vector& np_vec, const ProjData& stir) const; + + /// Convert STIR viewgram to SPECTGPU + void convert_viewgram_stir_to_SPECTGPU(std::vector& np_vec, const Viewgram& viewgram) const; + + /// Convert SPECTGPU proj data to STIR proj data + void convert_proj_data_SPECTGPU_to_stir(ProjData& stir_sptr, const std::vector& np_vec) const; + + /// Back project. Do some unavoidable const_casting as the wrapped methods don't use const + void back_project(std::vector& image, const std::vector& sino_no_gaps) const; + + /// Forward project, returns sinogram without gaps. Do some unavoidable const_casting as the wrapped methods don't use const + void forward_project(std::vector& sino_no_gaps, const std::vector& image) const; + + /// Create a STIR sinogram + static shared_ptr create_stir_sino(); + +private: + /// Check that set up has been run before returning data + void check_set_up() const; + + /// Permute the data + void permute(std::vector& output_array, + const std::vector& orig_array, + const unsigned output_dims[3], + const unsigned* permute_order) const; + + /// Convert 3d SPECTGPU proj data index to 1d + unsigned convert_SPECTGPU_proj_3d_to_1d_idx(const unsigned ang, const unsigned bins, const unsigned sino) const; + + /// Convert 1d SPECTGPU proj data index to 3d + void convert_SPECTGPU_proj_1d_to_3d_idx(unsigned& ang, unsigned& bins, unsigned& sino, const unsigned idx) const; + + bool _already_set_up; + char _devid; + shared_ptr _cnt_sptr; + int _nsinos; + char _att; + std::vector _isub; + bool _verbose; + Scanner::Type _scanner_type; + + + std::vector _crs; + std::vector _s2c; + +}; + +END_NAMESPACE_STIR + +#endif // __stir_gpu_SPECTGPUHelper_h__ diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index b0d6466df1..195c625120 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -135,6 +135,10 @@ if (STIR_WITH_CUDA) list(APPEND ${dir_LIB_SOURCES} CUDA/CudaRelativeDifferencePrior.cu CUDA/CudaGibbsPenaltiesInstantiations.cu + SPECTGPU_projector/BackProjectorByBinSPECTGPU.cxx + SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.cxx + SPECTGPU_projector/ProjectorByBinPairUsingSPECTGPU.cxx + SPECTGPU_projector/SPECTGPUHelper.cxx ) endif() diff --git a/src/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.cxx b/src/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.cxx new file mode 100644 index 0000000000..e02edc95f3 --- /dev/null +++ b/src/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.cxx @@ -0,0 +1,130 @@ +// +// +/*! + \file + \ingroup projection + \ingroup SPECTGPU + + \brief non-inline implementations for stir::BackProjectorByBinSPECTGPU + + \author Daniel Deidda + +*/ +/* + Copyright (C) 2026, National Physical Laboratory + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#include "stir/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.h" +#include "stir/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.h" +#include "stir/DiscretisedDensity.h" +#include "stir/RelatedViewgrams.h" +#include "stir/recon_buildblock/TrivialDataSymmetriesForBins.h" +#include "stir/ProjDataInfoCylindricalNoArcCorr.h" +#include "stir/recon_array_functions.h" +#include "stir/error.h" + +START_NAMESPACE_STIR + +////////////////////////////////////////////////////////// +const char* const BackProjectorByBinSPECTGPU::registered_name = "SPECTGPU"; + +BackProjectorByBinSPECTGPU::BackProjectorByBinSPECTGPU() + : _cuda_device(0), + _cuda_verbosity(true), + _use_truncation(false) +{ + this->_already_set_up = false; +} + +BackProjectorByBinSPECTGPU::~BackProjectorByBinSPECTGPU() +{} + +void +BackProjectorByBinSPECTGPU::initialise_keymap() +{ + parser.add_start_key("Back Projector Using SPECTGPU Parameters"); + parser.add_stop_key("End Back Projector Using SPECTGPU Parameters"); + parser.add_key("CUDA device", &_cuda_device); + parser.add_key("verbosity", &_cuda_verbosity); +} + +void +BackProjectorByBinSPECTGPU::set_up(const shared_ptr& proj_data_info_sptr, + const shared_ptr>& density_info_sptr) +{ + BackProjectorByBin::set_up(proj_data_info_sptr, density_info_sptr); + check(*proj_data_info_sptr, *_density_sptr); + _symmetries_sptr.reset(new TrivialDataSymmetriesForBins(proj_data_info_sptr)); + + + // Set up the SPECTGPU binary helper + _helper.set_cuda_device_id(_cuda_device); + _helper.set_scanner_type(proj_data_info_sptr->get_scanner_ptr()->get_type()); + _helper.set_att(0); + _helper.set_verbose(_cuda_verbosity); + _helper.set_up(); + + // Create sinogram + _np_sino = _helper.create_SPECTGPU_sinogram(); +} + +void +BackProjectorByBinSPECTGPU::back_project(const ProjData& proj_data, int subset_num, int num_subsets) +{ + // Check the user has tried to project all data + if (subset_num != 0 || num_subsets != 1) + error("BackProjectorByBinSPECTGPU::back_project " + "only works with all data (no subsets)."); + + _helper.convert_proj_data_stir_to_SPECTGPU(_np_sino, proj_data); +} + +void +BackProjectorByBinSPECTGPU::get_output(DiscretisedDensity<3, float>& density) const +{ + + std::vector sino = _helper.create_SPECTGPU_sinogram(); + + // --------------------------------------------------------------- // + // Back project + // --------------------------------------------------------------- // + + std::vector np_im = _helper.create_SPECTGPU_image(); + _helper.back_project(np_im, sino_); + + // --------------------------------------------------------------- // + // SPECTGPU -> STIR image conversion + // --------------------------------------------------------------- // + + _helper.convert_image_SPECTGPU_to_stir(density, np_im); + + // After the back projection, we enforce a truncation outside of the FOV. + // This is because the SPECTGPU FOV is smaller than the STIR FOV and this + // could cause some voxel values to spiral out of control. + if (_use_truncation) + truncate_rim(density, 17); +} + +void +BackProjectorByBinSPECTGPU::start_accumulating_in_new_target() +{ + // Call base level + BackProjectorByBin::start_accumulating_in_new_target(); + // Also reset the SPECTGPU sinogram + _np_sino = _helper.create_SPECTGPU_sinogram(); +} + +void +BackProjectorByBinSPECTGPU::actual_back_project( + const RelatedViewgrams& related_viewgrams, const int, const int, const int, const int) +{ + for (stir::RelatedViewgrams::const_iterator iter = related_viewgrams.begin(); iter != related_viewgrams.end(); ++iter) + _helper.convert_viewgram_stir_to_SPECTGPU(_np_sino_w_gaps, *iter); +} + +END_NAMESPACE_STIR diff --git a/src/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.cxx b/src/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.cxx new file mode 100644 index 0000000000..20431dacf9 --- /dev/null +++ b/src/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.cxx @@ -0,0 +1,132 @@ +// +// +/*! + + \file + \ingroup projection + \ingroup SPECTGPU + + \brief non-inline implementations for stir::ForwardProjectorByBinSPECTGPU + + \author Daniel Deidda + + +*/ +/* + Copyright (C) 2026, National Physical Laboratory + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#include "stir/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.h" +#include "stir/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.h" +#include "stir/ProjDataInMemory.h" +#include "stir/RelatedViewgrams.h" +#include "stir/ProjDataInfoCylindricalNoArcCorr.h" +#include "stir/recon_buildblock/TrivialDataSymmetriesForBins.h" +#include "stir/recon_array_functions.h" +#include "stir/error.h" + +START_NAMESPACE_STIR + +////////////////////////////////////////////////////////// +const char* const ForwardProjectorByBinSPECTGPU::registered_name = "SPECTGPU"; + +ForwardProjectorByBinSPECTGPU::ForwardProjectorByBinSPECTGPU() + : _cuda_device(0), + _cuda_verbosity(true), + _use_truncation(false) +{ + this->_already_set_up = false; +} + +ForwardProjectorByBinSPECTGPU::~ForwardProjectorByBinSPECTGPU() +{} + +void +ForwardProjectorByBinSPECTGPU::initialise_keymap() +{ + parser.add_start_key("Forward Projector Using SPECTGPU Parameters"); + parser.add_stop_key("End Forward Projector Using SPECTGPU Parameters"); + parser.add_key("CUDA device", &_cuda_device); + parser.add_key("verbosity", &_cuda_verbosity); +} + +void +ForwardProjectorByBinSPECTGPU::set_up(const shared_ptr& proj_data_info_sptr, + const shared_ptr>& density_info_sptr) +{ + ForwardProjectorByBin::set_up(proj_data_info_sptr, density_info_sptr); + check(*proj_data_info_sptr, *_density_sptr); + _symmetries_sptr.reset(new TrivialDataSymmetriesForBins(proj_data_info_sptr)); + + + // Initialise projected_data_sptr from this->_proj_data_info_sptr + _projected_data_sptr.reset(new ProjDataInMemory(this->_density_sptr->get_exam_info_sptr(), proj_data_info_sptr)); + + // Set up the SPECTGPU binary helper + _helper.set_scanner_type(proj_data_info_sptr->get_scanner_ptr()->get_type()); + _helper.set_cuda_device_id(_cuda_device); + _helper.set_att(0); + _helper.set_verbose(_cuda_verbosity); + _helper.set_up(); +} + + +void +ForwardProjectorByBinSPECTGPU::actual_forward_project( + RelatedViewgrams&, const DiscretisedDensity<3, float>&, const int, const int, const int, const int) +{ + throw std::runtime_error("Need to use set_input() if wanting to use ForwardProjectorByBinSPECTGPU."); +} + +void +ForwardProjectorByBinSPECTGPU::actual_forward_project( + RelatedViewgrams& viewgrams, const int, const int, const int, const int) +{ + // if (min_axial_pos_num != _proj_data_info_sptr->get_min_axial_pos_num() || + // ... ) + // error(); + + viewgrams = _projected_data_sptr->get_related_viewgrams(viewgrams.get_basic_view_segment_num(), _symmetries_sptr); +} + +void +ForwardProjectorByBinSPECTGPU::set_input(const DiscretisedDensity<3, float>& density) +{ + ForwardProjectorByBin::set_input(density); + + // Before forward projection, we enforce a truncation outside of the FOV. + // This is because the SPECTGPU FOV is smaller than the STIR FOV and this + // could cause some voxel values to spiral out of control. + if (_use_truncation) + truncate_rim(*_density_sptr, 17); + + // --------------------------------------------------------------- // + // STIR -> SPECTGPU image data conversion + // --------------------------------------------------------------- // + + std::vector np_vec = _helper.create_SPECTGPU_image(); + _helper.convert_image_stir_to_SPECTGPU(np_vec, *_density_sptr); + + // --------------------------------------------------------------- // + // Forward projection + // --------------------------------------------------------------- // + + std::vector sino = _helper.create_SPECTGPU_sinogram(); + _helper.forward_project(sino, np_vec); + + std::vector sino = _helper.create_SPECTGPU_sinogram(); + + + // --------------------------------------------------------------- // + // SPECTGPU -> STIR projection data conversion + // --------------------------------------------------------------- // + + _helper.convert_proj_data_SPECTGPU_to_stir(*_projected_data_sptr, sino); +} + +END_NAMESPACE_STIR diff --git a/src/recon_buildblock/SPECTGPU_projector/ProjectorByBinPairUsingSPECTGPU.cxx b/src/recon_buildblock/SPECTGPU_projector/ProjectorByBinPairUsingSPECTGPU.cxx new file mode 100644 index 0000000000..ad44e5ab2b --- /dev/null +++ b/src/recon_buildblock/SPECTGPU_projector/ProjectorByBinPairUsingSPECTGPU.cxx @@ -0,0 +1,98 @@ +// +// +/*! + \file + \ingroup projection + \ingroup SPECTGPU + + \brief non-inline implementations for stir::ProjectorByBinPairUsingSPECTGPU + + \author Daniel Deidda + +*/ +/* + Copyright (C) 2026, National Physical Laboratory + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#include "stir/recon_buildblock/SPECTGPU_projector/ProjectorByBinPairUsingSPECTGPU.h" +#include "stir/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.h" +#include "stir/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.h" + +START_NAMESPACE_STIR + +const char* const ProjectorByBinPairUsingSPECTGPU::registered_name = "SPECTGPU"; + +void +ProjectorByBinPairUsingSPECTGPU::initialise_keymap() +{ + base_type::initialise_keymap(); + parser.add_start_key("Projector Pair Using SPECTGPU Parameters"); + parser.add_stop_key("End Projector Pair Using SPECTGPU Parameters"); + parser.add_key("verbosity", &_verbosity); + parser.add_key("use_truncation", &_use_truncation); +} + +void +ProjectorByBinPairUsingSPECTGPU::set_defaults() +{ + base_type::set_defaults(); + this->set_verbosity(true); + this->set_use_truncation(false); +} + +bool +ProjectorByBinPairUsingSPECTGPU::post_processing() +{ + this->set_verbosity(this->_verbosity); + this->set_use_truncation(this->_use_truncation); + + if (base_type::post_processing()) + return true; + return false; +} + +ProjectorByBinPairUsingSPECTGPU::ProjectorByBinPairUsingSPECTGPU() +{ + this->forward_projector_sptr.reset(new ForwardProjectorByBinSPECTGPU); + this->back_projector_sptr.reset(new BackProjectorByBinSPECTGPU); + set_defaults(); +} + +void +ProjectorByBinPairUsingSPECTGPU::set_verbosity(const bool verbosity) +{ + _verbosity = verbosity; + + shared_ptr fwd_prj_downcast_sptr + = dynamic_pointer_cast(this->forward_projector_sptr); + if (fwd_prj_downcast_sptr) + fwd_prj_downcast_sptr->set_verbosity(_verbosity); + + shared_ptr bck_prj_downcast_sptr + = dynamic_pointer_cast(this->back_projector_sptr); + if (bck_prj_downcast_sptr) + bck_prj_downcast_sptr->set_verbosity(_verbosity); +} + +void +ProjectorByBinPairUsingSPECTGPU::set_use_truncation(const bool use_truncation) +{ + _use_truncation = use_truncation; + + shared_ptr fwd_prj_downcast_sptr + = dynamic_pointer_cast(this->forward_projector_sptr); + if (fwd_prj_downcast_sptr) + fwd_prj_downcast_sptr->set_use_truncation(_use_truncation); + + shared_ptr bck_prj_downcast_sptr + = dynamic_pointer_cast(this->back_projector_sptr); + if (bck_prj_downcast_sptr) + bck_prj_downcast_sptr->set_use_truncation(_use_truncation); +} + +END_NAMESPACE_STIR diff --git a/src/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.cxx b/src/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.cxx new file mode 100644 index 0000000000..2ca9788968 --- /dev/null +++ b/src/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.cxx @@ -0,0 +1,609 @@ +// +// +/*! + + \file + \ingroup projection + \ingroup SPECTGPU + + \brief non-inline implementations for stir::SPECTGPUHelper + + \author Daniel Deidda + + +*/ +/* + Copyright (C) 2026, National Physical Laboratory + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#include "stir/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.h" +#include "stir/VoxelsOnCartesianGrid.h" +#include "stir/is_null_ptr.h" +#include "stir/ProjDataInfoCylindricalNoArcCorr.h" +#include "stir/ProjDataInMemory.h" +#include "stir/IndexRange3D.h" +#include "stir/FilePath.h" +#include "stir/IO/stir_ecat_common.h" +#include "stir/error.h" +#include "stir/format.h" +// Non-STIR includes +#include +#include +#include "driver_types.h" +// SPECTGPU includes +#include "def.h" +#include "auxmath.h" +#include "prjb.h" +#include "prjf.h" +#include "recon.h" +#include "lmproc.h" +#include "scanner_0.h" +#include "rnd.h" +#include "norm.h" + +START_NAMESPACE_STIR + +SPECTGPUHelper::~SPECTGPUHelper() +{} + +// static void +// delete_axialLUT(axialLUT* axlut_ptr) +// { +// if (!axlut_ptr) +// return; +// delete[] axlut_ptr->li2rno; +// delete[] axlut_ptr->li2sn; +// delete[] axlut_ptr->li2nos; +// delete[] axlut_ptr->sn1_rno; +// delete[] axlut_ptr->sn1_sn11; +// delete[] axlut_ptr->sn1_ssrb; +// delete[] axlut_ptr->sn1_sn11no; +// } + +static shared_ptr +get_cnst(const Scanner& scanner, const bool cuda_verbose, const char cuda_device) +{ + shared_ptr cnt_sptr = MAKE_SHARED(); + + cnt_sptr->DEVID = cuda_device; // device (GPU) ID. allows choosing the device on which to perform calculations + cnt_sptr->VERBOSE = cuda_verbose; + + cnt_sptr->A = NSANGLES; // sino angles + cnt_sptr->W = NSBINS; // sino bins for any angular index + cnt_sptr->aw = AW; // sino bins (active only) + + cnt_sptr->NCRS = nCRS; // number of crystals + cnt_sptr->NRNG = NRINGS; // number of axial positions + cnt_sptr->D = -1; // number of linear indexes along Michelogram diagonals /*unknown*/ + cnt_sptr->Bt = -1; // number of buckets transaxially /*unknown*/ + + cnt_sptr->B = NBUCKTS; // number of buckets (total) + cnt_sptr->Cbt = 32552; // number of crystals in bucket transaxially /*unknown*/ + cnt_sptr->Cba = 3; // number of crystals in bucket axially /*unknown*/ + + cnt_sptr->NSN1 = NSINOS; // number of sinos + cnt_sptr->NSN64 = NRINGS * NRINGS; // with no MRD limit + cnt_sptr->NSEG0 = SEG0; + + cnt_sptr->RNG_STRT = 0; + cnt_sptr->RNG_END = NRINGS; + + cnt_sptr->ALPHA = aLPHA; // angle subtended by a crystal + float R = 32.8f; // ring radius + cnt_sptr->RE = R + 0.67f; // effective ring radius accounting for the depth of interaction + cnt_sptr->AXR = SZ_RING; // axial crystal dim + + float CLGHT = 29979245800.f; // speed of light [cm/s] + return cnt_sptr; +} + +static inline unsigned +to_1d_idx(const unsigned nrow, const unsigned ncol, const unsigned row, const unsigned col) +{ + return col + ncol * row; +} + +template +dataType* +create_heap_array(const unsigned numel, const dataType val = dataType(0)) +{ + dataType* array = new dataType[numel]; + std::fill(array, array + numel, val); + return array; +} + +void +SPECTGPUHelper::set_up() +{ + if (_att < 0) + throw std::runtime_error("SPECTGPUHelper::set_up() " + "emission or transmission mode (att) not set."); + + // Get consts + _cnt_sptr = get_cnst(_scanner_type, _verbose, _devid); + + + // isub + _isub = std::vector(unsigned(AW)); + for (unsigned i = 0; i < unsigned(AW); i++) + _isub[i] = int(i); + + _already_set_up = true; +} + +void +SPECTGPUHelper::check_set_up() const +{ + if (!_already_set_up) + throw std::runtime_error("SPECTGPUHelper::check_set_up() " + "Make sure filenames have been set and set_up has been run."); +} + +std::vector +SPECTGPUHelper::create_SPECTGPU_image() +{ + return std::vector(SZ_IMZ * SZ_IMX * SZ_IMY, 0); +} + +shared_ptr> +SPECTGPUHelper::create_stir_im() +{ + int nz(SZ_IMZ), nx(SZ_IMX), ny(SZ_IMY); + float sz(SZ_VOXZ * 10.f), sx(SZ_VOXY * 10.f), sy(SZ_VOXY * 10.f); + shared_ptr> out_im_stir_sptr = MAKE_SHARED>( + IndexRange3D(0, nz - 1, -(ny / 2), -(ny / 2) + ny - 1, -(nx / 2), -(nx / 2) + nx - 1), + CartesianCoordinate3D(0.f, 0.f, 0.f), + CartesianCoordinate3D(sz, sy, sx)); + return out_im_stir_sptr; +} + +std::vector +SPECTGPUHelper::create_SPECTGPU_sinogram() const +{ + check_set_up(); + return std::vector(_isub.size() * static_cast(_nsinos), 0); +} + + +void +get_stir_indices_and_dims(int stir_dim[3], + Coordinate3D& min_indices, + Coordinate3D& max_indices, + const DiscretisedDensity<3, float>& stir) +{ + if (!stir.get_regular_range(min_indices, max_indices)) + throw std::runtime_error("SPECTGPUHelper::set_input - " + "expected image to have regular range."); + for (int i = 0; i < 3; ++i) + stir_dim[i] = max_indices[i + 1] - min_indices[i + 1] + 1; +} + +unsigned +convert_SPECTGPU_im_3d_to_1d_idx(const unsigned x, const unsigned y, const unsigned z) +{ + return z * SZ_IMX * SZ_IMY + y * SZ_IMX + x; +} + +unsigned +SPECTGPUHelper::convert_SPECTGPU_proj_3d_to_1d_idx(const unsigned ang, const unsigned bins, const unsigned sino) const +{ + return sino * NSANGLES * NSBINS + ang * NSBINS + bins; +} + +void +SPECTGPUHelper::permute(std::vector& output_array, + const std::vector& orig_array, + const unsigned output_dims[3], + const unsigned permute_order[3]) const +{ +#ifndef NDEBUG + // Check that in the permute order, each number is between 0 and 2 (can't be <0 because it's unsigned) + for (unsigned i = 0; i < 3; ++i) + if (permute_order[i] > 2) + throw std::runtime_error("Permute order values should be between 0 and 2."); + // Check that each number is unique + for (unsigned i = 0; i < 3; ++i) + for (unsigned j = i + 1; j < 3; ++j) + if (permute_order[i] == permute_order[j]) + throw std::runtime_error("Permute order values should be unique."); + // Check that size of output_dims==arr.size() + assert(orig_array.size() == output_dims[0] * output_dims[1] * output_dims[2]); + // Check that output array is same size as input array + assert(orig_array.size() == output_array.size()); +#endif + + // Calculate old dimensions + unsigned old_dims[3]; + for (unsigned i = 0; i < 3; ++i) + old_dims[permute_order[i]] = output_dims[i]; + + // Loop over all elements + for (unsigned old_1d_idx = 0; old_1d_idx < orig_array.size(); ++old_1d_idx) + { + + // From the 1d index, generate the old 3d index + unsigned old_3d_idx[3] + = { old_1d_idx / (old_dims[2] * old_dims[1]), (old_1d_idx / old_dims[2]) % old_dims[1], old_1d_idx % old_dims[2] }; + + // Get the corresponding new 3d index + unsigned new_3d_idx[3]; + for (unsigned i = 0; i < 3; ++i) + new_3d_idx[i] = old_3d_idx[permute_order[i]]; + + // Get the new 1d index from the new 3d index + const unsigned new_1d_idx + = new_3d_idx[0] * output_dims[2] * output_dims[1] + new_3d_idx[1] * output_dims[2] + new_3d_idx[2]; + + // Fill the data + output_array[new_1d_idx] = orig_array[old_1d_idx]; + } +} + +void +SPECTGPUHelper::back_project(std::vector& image, const std::vector& sino_no_gaps) const +{ + check_set_up(); + assert(!image.empty()); + + std::vector unpermuted_image = this->create_SPECTGPU_image(); + + if (_verbose) + getMemUse(); + + // here start our code that calls the kernels + //rotation kernel + //backward kernel + +} + +void +SPECTGPUHelper::forward_project(std::vector& sino, const std::vector& image) const +{ + check_set_up(); + assert(!sino.empty()); + + // Permute the data (as this is done on the SPECTGPU python side before forward projection + // unsigned output_dims[3] = { 320, 320, 127 }; + // unsigned permute_order[3] = { 1, 2, 0 }; + // std::vector permuted_image = this->create_SPECTGPU_image(); + // this->permute(permuted_image, image, output_dims, permute_order); + + std::vector image = this->create_SPECTGPU_image(); + + if (_verbose) + getMemUse(); + + // here we do our calculation + // for view + // call rotation kernel + // call forward kernel +} + +shared_ptr +SPECTGPUHelper::create_stir_sino() +{ + + shared_ptr ei_sptr = MAKE_SHARED(); + ei_sptr->imaging_modality = ImagingModality::NM; + shared_ptr scanner_sptr(Scanner::get_scanner_from_name("uncknown")); + int num_views = scanner_sptr->get_num_detectors_per_ring() / 2 ; + int num_tang_pos = scanner_sptr->get_max_num_non_arccorrected_bins(); + shared_ptr pdi_sptr + = ProjDataInfo::construct_proj_data_info(scanner_sptr, span, max_ring_diff, num_views, num_tang_pos, false); + shared_ptr pd_sptr = MAKE_SHARED(ei_sptr, pdi_sptr); + return pd_sptr; +} + +template +static dataType* +read_from_binary_file(std::ifstream& file, const unsigned long num_elements) +{ + // Get current position, get size to end and go back to current position + const unsigned long current_pos = file.tellg(); + file.seekg(std::ios::cur, std::ios::end); + const unsigned long remaining_elements = file.tellg() / sizeof(dataType); + file.seekg(current_pos, std::ios::beg); + + if (remaining_elements < num_elements) + throw std::runtime_error("File smaller than requested."); + + dataType* contents = create_heap_array(num_elements); + file.read(reinterpret_cast(contents), num_elements * sizeof(dataType)); + return contents; +} + + +void +check_im_sizes(const int stir_dim[3], const int np_dim[3]) +{ + for (int i = 0; i < 3; ++i) + if (stir_dim[i] != np_dim[i]) + throw std::runtime_error(format("SPECTGPUHelper::check_im_sizes() - " + "STIR image ({}, {}, {}) should be == ({},{},{}).", + stir_dim[0], + stir_dim[1], + stir_dim[2], + np_dim[0], + np_dim[1], + np_dim[2])); +} + +void +check_voxel_spacing(const DiscretisedDensity<3, float>& stir) +{ + // Requires image to be a VoxelsOnCartesianGrid + const VoxelsOnCartesianGrid& stir_vocg = dynamic_cast&>(stir); + const BasicCoordinate<3, float> stir_spacing = stir_vocg.get_grid_spacing(); + + // Get SPECTGPU image spacing (need to *10 for mm) + float np_spacing[3] = { 10.f * SZ_VOXZ, 10.f * SZ_VOXY, 10.f * SZ_VOXY }; + + for (unsigned i = 0; i < 3; ++i) + if (std::abs(stir_spacing[int(i) + 1] - np_spacing[i]) > 1e-4f) + throw std::runtime_error(format("SPECTGPUHelper::check_voxel_spacing() - " + "STIR image ({}, {}, {}) should be == ({},{},{}).", + stir_spacing[1], + stir_spacing[2], + stir_spacing[3], + np_spacing[0], + np_spacing[1], + np_spacing[2])); +} + +void +SPECTGPUHelper::convert_image_stir_to_SPECTGPU(std::vector& np_vec, const DiscretisedDensity<3, float>& stir) +{ + // Get the dimensions of the input image + Coordinate3D min_indices; + Coordinate3D max_indices; + int stir_dim[3]; + get_stir_indices_and_dims(stir_dim, min_indices, max_indices, stir); + + // SPECTGPU requires the image to be (z,x,y)=(SZ_IMZ,SZ_IMX,SZ_IMY) + // which at the time of writing was (127,320,320). + const int np_dim[3] = { SZ_IMZ, SZ_IMX, SZ_IMY }; + check_im_sizes(stir_dim, np_dim); + check_voxel_spacing(stir); + + // Copy data from STIR to SPECTGPU image + unsigned np_z, np_y, np_x, np_1d; + for (int z = min_indices[1]; z <= max_indices[1]; z++) + { + for (int y = min_indices[2]; y <= max_indices[2]; y++) + { + for (int x = min_indices[3]; x <= max_indices[3]; x++) + { + // Convert the stir 3d index to a SPECTGPU 1d index + np_z = unsigned(z - min_indices[1]); + np_y = unsigned(y - min_indices[2]); + np_x = unsigned(x - min_indices[3]); + np_1d = convert_SPECTGPU_im_3d_to_1d_idx(np_x, np_y, np_z); + np_vec[np_1d] = stir[z][y][x]; + } + } + } +} + +void +SPECTGPUHelper::convert_image_SPECTGPU_to_stir(DiscretisedDensity<3, float>& stir, const std::vector& np_vec) +{ + // Get the dimensions of the input image + Coordinate3D min_indices; + Coordinate3D max_indices; + int stir_dim[3]; + get_stir_indices_and_dims(stir_dim, min_indices, max_indices, stir); + + // SPECTGPU requires the image to be (z,x,y)=(SZ_IMZ,SZ_IMX,SZ_IMY) + // which at the time of writing was (127,320,320). + const int np_dim[3] = { SZ_IMZ, SZ_IMX, SZ_IMY }; + check_im_sizes(stir_dim, np_dim); + check_voxel_spacing(stir); + + // Copy data from SPECTGPU to STIR image + unsigned np_z, np_y, np_x, np_1d; + for (int z = min_indices[1]; z <= max_indices[1]; z++) + { + for (int y = min_indices[2]; y <= max_indices[2]; y++) + { + for (int x = min_indices[3]; x <= max_indices[3]; x++) + { + // Convert the stir 3d index to a SPECTGPU 1d index + np_z = unsigned(z - min_indices[1]); + np_y = unsigned(y - min_indices[2]); + np_x = unsigned(x - min_indices[3]); + np_1d = convert_SPECTGPU_im_3d_to_1d_idx(np_x, np_y, np_z); + stir[z][y][x] = np_vec[np_1d]; + } + } + } +} + +void +get_vals_for_proj_data_conversion(std::vector& sizes, + std::vector& segment_sequence, + int& num_sinograms, + int& min_view, + int& max_view, + int& min_tang_pos, + int& max_tang_pos, + const ProjDataInfo& proj_data_info, + const std::vector& np_vec) +{ + const ProjDataInfoCylindricalNoArcCorr* info_sptr = dynamic_cast(&proj_data_info); + if (is_null_ptr(info_sptr)) + error("SPECTGPUHelper: only works with cylindrical projection data without arc-correction"); + + segment_sequence = ecat::find_segment_sequence(proj_data_info); + sizes.resize(segment_sequence.size()); + for (std::size_t s = 0U; s < segment_sequence.size(); ++s) + sizes[s] = proj_data_info.get_num_axial_poss(segment_sequence[s]); + + // Get dimensions of STIR sinogram + min_view = proj_data_info.get_min_view_num(); + max_view = proj_data_info.get_max_view_num(); + min_tang_pos = proj_data_info.get_min_tangential_pos_num(); + max_tang_pos = proj_data_info.get_max_tangential_pos_num(); + + num_sinograms = proj_data_info.get_num_axial_poss(0); + for (int s = 1; s <= proj_data_info.get_max_segment_num(); ++s) + num_sinograms += 2 * proj_data_info.get_num_axial_poss(s); + + int num_proj_data_elems = num_sinograms * (1 + max_view - min_view) * (1 + max_tang_pos - min_tang_pos); + + // Make sure they're the same size + if (np_vec.size() != unsigned(num_proj_data_elems)) + error(format("SPECTGPUHelper::get_vals_for_proj_data_conversion " + "SPECTGPU and STIR sinograms are different sizes ({} for STIR versus {} for NP", + num_proj_data_elems, + np_vec.size())); +} + +void +get_stir_segment_and_axial_pos_from_SPECTGPU_sino( + int& segment, int& axial_pos, const unsigned np_sino, const std::vector& sizes, const std::vector& segment_sequence) +{ + int z = int(np_sino); + for (unsigned i = 0; i < segment_sequence.size(); ++i) + { + if (z < sizes[i]) + { + axial_pos = z; + segment = segment_sequence[i]; + return; + } + else + { + z -= sizes[i]; + } + } +} + +void +get_SPECTGPU_sino_from_stir_segment_and_axial_pos(unsigned& np_sino, + const int segment, + const int axial_pos, + const std::vector& sizes, + const std::vector& segment_sequence) +{ + np_sino = 0U; + for (unsigned i = 0; i < segment_sequence.size(); ++i) + { + if (segment == segment_sequence[i]) + { + np_sino += axial_pos; + return; + } + else + { + np_sino += sizes[i]; + } + } + throw std::runtime_error( + "SPECTGPUHelper::get_SPECTGPU_sino_from_stir_segment_and_axial_pos(): Failed to find SPECTGPU sinogram."); +} + +void +SPECTGPUHelper::convert_viewgram_stir_to_SPECTGPU(std::vector& np_vec, const Viewgram& viewgram) const +{ + // Get the values (and LUT) to be able to switch between STIR and SPECTGPU projDatas + std::vector sizes, segment_sequence; + int num_sinograms, min_view, max_view, min_tang_pos, max_tang_pos; + get_vals_for_proj_data_conversion(sizes, + segment_sequence, + num_sinograms, + min_view, + max_view, + min_tang_pos, + max_tang_pos, + *viewgram.get_proj_data_info_sptr(), + np_vec); + + const int segment = viewgram.get_segment_num(); + const int view = viewgram.get_view_num(); + + // Loop over the STIR view and tangential position + for (int ax_pos = viewgram.get_min_axial_pos_num(); ax_pos <= viewgram.get_max_axial_pos_num(); ++ax_pos) + { + + unsigned np_sino; + + // Convert the SPECTGPU sinogram to STIR's segment and axial position + get_SPECTGPU_sino_from_stir_segment_and_axial_pos(np_sino, segment, ax_pos, sizes, segment_sequence); + + for (int tang_pos = min_tang_pos; tang_pos <= max_tang_pos; ++tang_pos) + { + + unsigned np_ang = unsigned(view - min_view); + unsigned np_bin = unsigned(tang_pos - min_tang_pos); + unsigned np_1d = convert_SPECTGPU_proj_3d_to_1d_idx(np_ang, np_bin, np_sino); + np_vec.at(np_1d) = viewgram.at(ax_pos).at(tang_pos); + } + } +} + +void +SPECTGPUHelper::convert_proj_data_stir_to_SPECTGPU(std::vector& np_vec, const ProjData& stir) const +{ + const int min_view = stir.get_min_view_num(); + const int max_view = stir.get_max_view_num(); + const int min_segment = stir.get_min_segment_num(); + const int max_segment = stir.get_max_segment_num(); + + for (int view = min_view; view <= max_view; ++view) + { + for (int segment = min_segment; segment <= max_segment; ++segment) + { + convert_viewgram_stir_to_SPECTGPU(np_vec, stir.get_viewgram(view, segment)); + } + } +} + +void +SPECTGPUHelper::convert_proj_data_SPECTGPU_to_stir(ProjData& stir, const std::vector& np_vec) const +{ + // Get the values (and LUT) to be able to switch between STIR and SPECTGPU projDatas + std::vector sizes, segment_sequence; + int num_sinograms, min_view, max_view, min_tang_pos, max_tang_pos; + get_vals_for_proj_data_conversion(sizes, + segment_sequence, + num_sinograms, + min_view, + max_view, + min_tang_pos, + max_tang_pos, + *stir.get_proj_data_info_sptr(), + np_vec); + + int segment, axial_pos; + // Loop over all SPECTGPU sinograms + for (unsigned np_sino = 0; np_sino < unsigned(num_sinograms); ++np_sino) + { + + // Convert the SPECTGPU sinogram to STIR's segment and axial position + get_stir_segment_and_axial_pos_from_SPECTGPU_sino(segment, axial_pos, np_sino, sizes, segment_sequence); + + // Get the corresponding STIR sinogram + Sinogram sino = stir.get_empty_sinogram(axial_pos, segment); + + // Loop over the STIR view and tangential position + for (int view = min_view; view <= max_view; ++view) + { + for (int tang_pos = min_tang_pos; tang_pos <= max_tang_pos; ++tang_pos) + { + + unsigned np_ang = unsigned(view - min_view); + unsigned np_bin = unsigned(tang_pos - min_tang_pos); + unsigned np_1d = convert_SPECTGPU_proj_3d_to_1d_idx(np_ang, np_bin, np_sino); + sino.at(view).at(tang_pos) = np_vec.at(np_1d); + } + } + stir.set_sinogram(sino); + } +} + +END_NAMESPACE_STIR diff --git a/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu b/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu new file mode 100644 index 0000000000..8ea094d6e8 --- /dev/null +++ b/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu @@ -0,0 +1,52 @@ +// +// +/*! + + \file + \ingroup projection + \ingroup SPECTGPU + + \brief implementations for cuda kernel for rotating projector with gaussian interpolation + a la Wallis et al 1997,TMI, doi: 10.1109/42.552061. + + \author Daniel Deidda + + +*/ +/* + Copyright (C) 2026, National Physical Laboratory + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#include "stir/cuda_utilities.h" +#include +#include + + +//the following is a pull operation +__global__ void forwardKernel(const float* __restrict__ in_im, + float* __restrict__ out_sino, + int3 dim, + float3 spacing, + float3 origin, + float angle_rad) + { + + + } + +// the following is the adjoint operation (push) +__global__ void backwardKernel(const float* __restrict__ in_sino, + float* __restrict__ out_im, + int3 dim, + float3 spacing, + float3 origin, + float angle_rad) + { + + + } diff --git a/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu b/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu new file mode 100644 index 0000000000..2f2bc7830c --- /dev/null +++ b/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu @@ -0,0 +1,51 @@ +// +// +/*! + + \file + \ingroup projection + \ingroup SPECTGPU + + \brief implementations for cuda kernel for rotating projector with gaussian interpolation + a la Wallis et al 1997,TMI, doi: 10.1109/42.552061. + + \author Daniel Deidda + + +*/ +/* + Copyright (C) 2026, National Physical Laboratory + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#include "stir/cuda_utilities.h" +#include +#include + +//the following is a pull operation +__global__ void rotateKernel_pull(const float* __restrict__ in_im, + float* __restrict__ out_im, + int3 dim, + float3 spacing, + float3 origin, + float angle_rad) + { + + + } + +// the following is the adjoint operation (push) +__global__ void rotateKernel_push(const float* __restrict__ in_im, + float* __restrict__ out_im, + int3 dim, + float3 spacing, + float3 origin, + float angle_rad) + { + + + } From c948b16b48c48214b19d98e7f5317bbfe444ca25 Mon Sep 17 00:00:00 2001 From: williamuic <1187541557@qq.com> Date: Fri, 6 Mar 2026 14:13:58 +0000 Subject: [PATCH 2/9] forward line integrals --- .../SPECTGPU_projector/SPECTGPUProjection.cu | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu b/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu index 8ea094d6e8..f2a82628db 100644 --- a/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu +++ b/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu @@ -35,10 +35,22 @@ __global__ void forwardKernel(const float* __restrict__ in_im, float3 origin, float angle_rad) { + //1. define position in the detector, x horizontal, z vertical, y is the depth. + + int det_x = blockIdx.x * blockDim.x + threadIdx.x; + int det_y = blockIdx.y * blockDim.y + threadIdx.y; + int det_z = blockIdx.z * blockDim.z + threadIdx.z; + + if (det_x >= dim.x || det_y >= dim.y || det_z >= dim.z) { + return; + } + int sino_idx = det_z * dim.x + det_x; - - } - + // Sum voxel values along the y-axis (depth) for this detector pixel + int voxel_idx = det_y * dim.x * dim.z + det_z * dim.x + det_x; + out_sino[sino_idx] += in_im[voxel_idx]; + + } // the following is the adjoint operation (push) __global__ void backwardKernel(const float* __restrict__ in_sino, float* __restrict__ out_im, From 5d300cf8d59ed9a064370e04cb8f3d2827ea7f4f Mon Sep 17 00:00:00 2001 From: "JCHAN\\jchan" Date: Fri, 6 Mar 2026 10:45:57 -0500 Subject: [PATCH 3/9] implemented a backward projection kernel in cuda --- .../SPECTGPU_projector/SPECTGPUProjection.cu | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu b/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu index f2a82628db..6aa040140b 100644 --- a/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu +++ b/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu @@ -59,6 +59,11 @@ __global__ void backwardKernel(const float* __restrict__ in_sino, float3 origin, float angle_rad) { - - - } + int idz = blockIdx.zblockDim.z + threadIdx.z; + int idy = blockIdx.yblockDim.y + threadIdx.y; + int idx = blockIdx.xblockDim.x + threadIdx.x; + if ( idz < dim.z && idy < dim.y && idx < dim.x) + { + out_im[(idz dim.y + idy) * dim.x + idx] += sinogram[idz * dim.y + idy] * spacing.x; + } +} From d4cd638c8a3ceac86f1b9507fb12cbbff44eb5ee Mon Sep 17 00:00:00 2001 From: "david.roddy" Date: Fri, 6 Mar 2026 15:48:24 +0000 Subject: [PATCH 4/9] GPUSPECT: forward rotation and interpolation in cuda --- .../SPECTGPURotateAndGaussianInterpolate.cu | 99 ++++++++++++++++++- 1 file changed, 98 insertions(+), 1 deletion(-) diff --git a/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu b/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu index 2f2bc7830c..7bfc9f2b5a 100644 --- a/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu +++ b/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu @@ -34,9 +34,106 @@ __global__ void rotateKernel_pull(const float* __restrict__ in_im, float3 origin, float angle_rad) { + // parallelise the operation across all image voxels + int i = threadIdx.x + blockDim.x * blockIdx.x; + int j = threadIdx.y + blockDim.y * blockIdx.y; + int k = threadIdx.z + blockDim.z * blockIdx.z; + // check we are not outside the image + if (i >= dim.x || j >= dim.y || k >= dim.z) { + return; + }; - } + // get the 1-dimensional index + int idx = i + (j * dim.x) + (k * dim.x * dim.y); + + // move to centre of voxel and convert from voxel to image space + float Xs_rel_x = (i + 0.5f) * spacing.x - origin.x; + float Xs_rel_y = (j + 0.5f) * spacing.y - origin.y; + float Xs_rel_z = (k + 0.5f) * spacing.z - origin.z; + + // perform the rotation by 'angle_rad' + float X_rot_x = (Xs_rel_x * cosf(angle_rad)) - (Xs_rel_y * sinf(angle_rad)); + float X_rot_y = (Xs_rel_x * sinf(angle_rad)) + (Xs_rel_y * cosf(angle_rad)); + float X_rot_z = Xs_rel_z; + + // convert back to voxel space origin (corner of image) but keep image space dimensions + float Xr_x = X_rot_x + origin.x; + float Xr_y = X_rot_y + origin.y; + float Xr_z = X_rot_z + origin.z; + + // now we go fully back to voxel space (but now it's been rotated) + float p_r = Xr_x / spacing.x; + float q_r = Xr_y / spacing.y; + float r_r = Xr_z / spacing.z; + + // first loop is for finding the value of a gaussian kernel at + // each nearest neighbour position and summing them all. This will allow normalisation + // of each NN contribution in the next loop. + float G = 0; // accumulator variable + float sigma = 1; // gaussian kernel sigma + + for (int dr = -1; dr <= 1; dr++) { + int r = (int)roundf(r_r) + dr; // nearest neighbour z coordinate in rotated voxel space + float x_r = (r * spacing.z) - origin.z; // converted to image space + + for (int dq = -1; dq <= 1; dq++) { + int q = (int)roundf(q_r) + dq; + float x_q = (q * spacing.y) - origin.y; + + for (int dp = -1; dp <= 1; dp++) { + int p = (int)roundf(p_r) + dp; + float x_p = (p * spacing.x) - origin.x; + + // get distance between nearest neighbour and central voxel + // both need to be in image space + float delta_x = X_rot_x - x_p; + float delta_y = X_rot_y - x_q; + float delta_z = X_rot_z - x_r; + + // caulculate gaussian kernel + float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma)) + G += g; + }; + }; + }; + + // loop again but this time actually fetch the image values + float accumulation = 0; + + for (int dr = -1; dr <= 1; dr++) { + float r = dr + (int)roundf(r_r); + float x_r = origin.z + (r * spacing.z); + + for (int dq = -1; dq <= 1; dq++) { + float q = dq + (int)roundf(q_r); + float x_q = origin.y + (q * spacing.y); + + for (int dp = -1; dp <= 1; dp++) { + float p = dp + (int)roundf(p_r); + float x_p = origin.x + (p * spacing.x); + + float delta_x = Xr_x - x_p; + float delta_y = Xr_y - x_q; + float delta_z = Xr_z - x_r; + + // get input image voxel index that we are at + int i_idx = p + (q * dim.x) + (r * dim.x * dim.y); + + // calculate gaussian weighting for this voxel + float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma)); + + // record the weighted value from this voxel + accumulation += in_im[i_idx] * g / G; + + }; + }; + }; + + // assign the pulled voxel values to a single voxel in the output image + out_im[idx] = accumulation; + +}; // the following is the adjoint operation (push) __global__ void rotateKernel_push(const float* __restrict__ in_im, From 9ca0e11ba9c2e45f7541ce841b06d3f11e357f28 Mon Sep 17 00:00:00 2001 From: Catherine Buck Date: Fri, 6 Mar 2026 15:58:09 +0000 Subject: [PATCH 5/9] Adding kernel push --- .../SPECTGPURotateAndGaussianInterpolate.cu | 95 ++++++++++++++++++- 1 file changed, 91 insertions(+), 4 deletions(-) diff --git a/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu b/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu index 7bfc9f2b5a..5f168ccf90 100644 --- a/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu +++ b/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu @@ -40,9 +40,7 @@ __global__ void rotateKernel_pull(const float* __restrict__ in_im, int k = threadIdx.z + blockDim.z * blockIdx.z; // check we are not outside the image - if (i >= dim.x || j >= dim.y || k >= dim.z) { - return; - }; + if (i >= dim.x || j >= dim.y || k >= dim.z) return; // get the 1-dimensional index int idx = i + (j * dim.x) + (k * dim.x * dim.y); @@ -143,6 +141,95 @@ __global__ void rotateKernel_push(const float* __restrict__ in_im, float3 origin, float angle_rad) { + // parallelise the operation across all image voxels + int i = threadIdx.x + blockDim.x * blockIdx.x; + int j = threadIdx.y + blockDim.y * blockIdx.y; + int k = threadIdx.z + blockDim.z * blockIdx.z; + + // check we are not outside the image + if (i >= dim.x || j >= dim.y || k >= dim.z) return; + + // get the 1-dimensional index + int idx = i + (j * dim.x) + (k * dim.x * dim.y); + + // move to centre of voxel and convert from voxel to image space + float Xs_rel_x = (i + 0.5f) * spacing.x - origin.x; + float Xs_rel_y = (j + 0.5f) * spacing.y - origin.y; + float Xs_rel_z = (k + 0.5f) * spacing.z - origin.z; + + // perform the rotation by 'angle_rad' + float X_rot_x = (Xs_rel_x * cosf(angle_rad)) - (Xs_rel_y * sinf(angle_rad)); + float X_rot_y = (Xs_rel_x * sinf(angle_rad)) + (Xs_rel_y * cosf(angle_rad)); + float X_rot_z = Xs_rel_z; + + // convert back to voxel space origin (corner of image) but keep image space dimensions + float Xr_x = X_rot_x + origin.x; + float Xr_y = X_rot_y + origin.y; + float Xr_z = X_rot_z + origin.z; + + // now we go fully back to voxel space (but now it's been rotated) + float p_r = Xr_x / spacing.x; + float q_r = Xr_y / spacing.y; + float r_r = Xr_z / spacing.z; + + // first loop is for finding the value of a gaussian kernel at + // each nearest neighbour position and summing them all. This will allow normalisation + // of each NN contribution in the next loop. + float G = 0; // accumulator variable + float sigma = 1; // gaussian kernel sigma + + for (int dr = -1; dr <= 1; dr++) { + int r = (int)roundf(r_r) + dr; // nearest neighbour z coordinate in rotated voxel space + float x_r = (r * spacing.z) - origin.z; // converted to image space + + for (int dq = -1; dq <= 1; dq++) { + int q = (int)roundf(q_r) + dq; + float x_q = (q * spacing.y) - origin.y; + + for (int dp = -1; dp <= 1; dp++) { + int p = (int)roundf(p_r) + dp; + float x_p = (p * spacing.x) - origin.x; + + // get distance between nearest neighbour and central voxel + // both need to be in image space + float delta_x = X_rot_x - x_p; + float delta_y = X_rot_y - x_q; + float delta_z = X_rot_z - x_r; + + // caulculate gaussian kernel + float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma)) + G += g; + }; + }; + }; + + // loop again but this time actually fetch the image values + float accumulation = 0; + for (int dr = -1; dr <= 1; dr++) { + float r = dr + (int)roundf(r_r); + float x_r = origin.z + (r * spacing.z); + + for (int dq = -1; dq <= 1; dq++) { + float q = dq + (int)roundf(q_r); + float x_q = origin.y + (q * spacing.y); + + for (int dp = -1; dp <= 1; dp++) { + float p = dp + (int)roundf(p_r); + float x_p = origin.x + (p * spacing.x); - } + float delta_x = Xr_x - x_p; + float delta_y = Xr_y - x_q; + float delta_z = Xr_z - x_r; + + // get input image voxel index that we are at + int i_idx = p + (q * dim.x) + (r * dim.x * dim.y); + + // Pushing the weighted counts to NN + float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma)); + atomic_add(&out_im[i_idx], in_im[idx]*g/G); + + }; + }; + }; +} From d731b5c5728c9d7680d07629d5e7dd08696a21e3 Mon Sep 17 00:00:00 2001 From: Catherine Buck Date: Fri, 6 Mar 2026 16:04:13 +0000 Subject: [PATCH 6/9] fixup! Adding kernel push --- .../SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu b/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu index 5f168ccf90..2b6bff150e 100644 --- a/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu +++ b/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu @@ -227,7 +227,7 @@ __global__ void rotateKernel_push(const float* __restrict__ in_im, // Pushing the weighted counts to NN float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma)); - atomic_add(&out_im[i_idx], in_im[idx]*g/G); + atomicAdd(&out_im[i_idx], in_im[idx]*g/G); }; }; From 8dfbdd11c290ad3bc823c91648ef469f908a9b71 Mon Sep 17 00:00:00 2001 From: "david.roddy" Date: Fri, 6 Mar 2026 16:30:43 +0000 Subject: [PATCH 7/9] SPECTGPU: run pre-commit hooks on SPECTGPURotateAndGaussianInterpolate --- .../SPECTGPURotateAndGaussianInterpolate.cu | 241 +++++++++--------- 1 file changed, 123 insertions(+), 118 deletions(-) diff --git a/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu b/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu index 2b6bff150e..1e19382d67 100644 --- a/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu +++ b/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu @@ -6,7 +6,7 @@ \ingroup projection \ingroup SPECTGPU - \brief implementations for cuda kernel for rotating projector with gaussian interpolation + \brief implementations for cuda kernel for rotating projector with gaussian interpolation a la Wallis et al 1997,TMI, doi: 10.1109/42.552061. \author Daniel Deidda @@ -26,21 +26,19 @@ #include #include -//the following is a pull operation -__global__ void rotateKernel_pull(const float* __restrict__ in_im, - float* __restrict__ out_im, - int3 dim, - float3 spacing, - float3 origin, - float angle_rad) - { +// the following is a pull operation +__global__ void +rotateKernel_pull( + const float* __restrict__ in_im, float* __restrict__ out_im, int3 dim, float3 spacing, float3 origin, float angle_rad) +{ // parallelise the operation across all image voxels int i = threadIdx.x + blockDim.x * blockIdx.x; int j = threadIdx.y + blockDim.y * blockIdx.y; int k = threadIdx.z + blockDim.z * blockIdx.z; // check we are not outside the image - if (i >= dim.x || j >= dim.y || k >= dim.z) return; + if (i >= dim.x || j >= dim.y || k >= dim.z) + return; // get the 1-dimensional index int idx = i + (j * dim.x) + (k * dim.x * dim.y); @@ -68,86 +66,88 @@ __global__ void rotateKernel_pull(const float* __restrict__ in_im, // first loop is for finding the value of a gaussian kernel at // each nearest neighbour position and summing them all. This will allow normalisation // of each NN contribution in the next loop. - float G = 0; // accumulator variable - float sigma = 1; // gaussian kernel sigma - - for (int dr = -1; dr <= 1; dr++) { - int r = (int)roundf(r_r) + dr; // nearest neighbour z coordinate in rotated voxel space - float x_r = (r * spacing.z) - origin.z; // converted to image space - - for (int dq = -1; dq <= 1; dq++) { - int q = (int)roundf(q_r) + dq; - float x_q = (q * spacing.y) - origin.y; - - for (int dp = -1; dp <= 1; dp++) { - int p = (int)roundf(p_r) + dp; - float x_p = (p * spacing.x) - origin.x; - - // get distance between nearest neighbour and central voxel - // both need to be in image space - float delta_x = X_rot_x - x_p; - float delta_y = X_rot_y - x_q; - float delta_z = X_rot_z - x_r; - - // caulculate gaussian kernel - float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma)) - G += g; - }; + float G = 0; // accumulator variable + float sigma = 1; // gaussian kernel sigma + + for (int dr = -1; dr <= 1; dr++) + { + int r = (int)roundf(r_r) + dr; // nearest neighbour z coordinate in rotated voxel space + float x_r = (r * spacing.z) - origin.z; // converted to image space + + for (int dq = -1; dq <= 1; dq++) + { + int q = (int)roundf(q_r) + dq; + float x_q = (q * spacing.y) - origin.y; + + for (int dp = -1; dp <= 1; dp++) + { + int p = (int)roundf(p_r) + dp; + float x_p = (p * spacing.x) - origin.x; + + // get distance between nearest neighbour and central voxel + // both need to be in image space + float delta_x = X_rot_x - x_p; + float delta_y = X_rot_y - x_q; + float delta_z = X_rot_z - x_r; + + // caulculate gaussian kernel + float g = expf(-1. * ((delta_x * delta_x) + (delta_y * delta_y) + (delta_z * delta_z)) / (2 * sigma * sigma)) G + += g; + }; + }; }; - }; // loop again but this time actually fetch the image values float accumulation = 0; - for (int dr = -1; dr <= 1; dr++) { - float r = dr + (int)roundf(r_r); - float x_r = origin.z + (r * spacing.z); + for (int dr = -1; dr <= 1; dr++) + { + float r = dr + (int)roundf(r_r); + float x_r = origin.z + (r * spacing.z); - for (int dq = -1; dq <= 1; dq++) { - float q = dq + (int)roundf(q_r); - float x_q = origin.y + (q * spacing.y); + for (int dq = -1; dq <= 1; dq++) + { + float q = dq + (int)roundf(q_r); + float x_q = origin.y + (q * spacing.y); - for (int dp = -1; dp <= 1; dp++) { - float p = dp + (int)roundf(p_r); - float x_p = origin.x + (p * spacing.x); + for (int dp = -1; dp <= 1; dp++) + { + float p = dp + (int)roundf(p_r); + float x_p = origin.x + (p * spacing.x); - float delta_x = Xr_x - x_p; - float delta_y = Xr_y - x_q; - float delta_z = Xr_z - x_r; + float delta_x = Xr_x - x_p; + float delta_y = Xr_y - x_q; + float delta_z = Xr_z - x_r; - // get input image voxel index that we are at - int i_idx = p + (q * dim.x) + (r * dim.x * dim.y); + // get input image voxel index that we are at + int i_idx = p + (q * dim.x) + (r * dim.x * dim.y); - // calculate gaussian weighting for this voxel - float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma)); + // calculate gaussian weighting for this voxel + float g = expf(-1. * ((delta_x * delta_x) + (delta_y * delta_y) + (delta_z * delta_z)) / (2 * sigma * sigma)); - // record the weighted value from this voxel - accumulation += in_im[i_idx] * g / G; - - }; + // record the weighted value from this voxel + accumulation += in_im[i_idx] * g / G; + }; + }; }; - }; // assign the pulled voxel values to a single voxel in the output image out_im[idx] = accumulation; - }; // the following is the adjoint operation (push) -__global__ void rotateKernel_push(const float* __restrict__ in_im, - float* __restrict__ out_im, - int3 dim, - float3 spacing, - float3 origin, - float angle_rad) - { - // parallelise the operation across all image voxels +__global__ void +rotateKernel_push( + const float* __restrict__ in_im, float* __restrict__ out_im, int3 dim, float3 spacing, float3 origin, float angle_rad) +{ + // parallelise the operation across all image voxels int i = threadIdx.x + blockDim.x * blockIdx.x; int j = threadIdx.y + blockDim.y * blockIdx.y; int k = threadIdx.z + blockDim.z * blockIdx.z; // check we are not outside the image - if (i >= dim.x || j >= dim.y || k >= dim.z) return; + if (i >= dim.x || j >= dim.y || k >= dim.z) + return; // get the 1-dimensional index int idx = i + (j * dim.x) + (k * dim.x * dim.y); @@ -175,61 +175,66 @@ __global__ void rotateKernel_push(const float* __restrict__ in_im, // first loop is for finding the value of a gaussian kernel at // each nearest neighbour position and summing them all. This will allow normalisation // of each NN contribution in the next loop. - float G = 0; // accumulator variable - float sigma = 1; // gaussian kernel sigma - - for (int dr = -1; dr <= 1; dr++) { - int r = (int)roundf(r_r) + dr; // nearest neighbour z coordinate in rotated voxel space - float x_r = (r * spacing.z) - origin.z; // converted to image space - - for (int dq = -1; dq <= 1; dq++) { - int q = (int)roundf(q_r) + dq; - float x_q = (q * spacing.y) - origin.y; - - for (int dp = -1; dp <= 1; dp++) { - int p = (int)roundf(p_r) + dp; - float x_p = (p * spacing.x) - origin.x; - - // get distance between nearest neighbour and central voxel - // both need to be in image space - float delta_x = X_rot_x - x_p; - float delta_y = X_rot_y - x_q; - float delta_z = X_rot_z - x_r; - - // caulculate gaussian kernel - float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma)) - G += g; - }; + float G = 0; // accumulator variable + float sigma = 1; // gaussian kernel sigma + + for (int dr = -1; dr <= 1; dr++) + { + int r = (int)roundf(r_r) + dr; // nearest neighbour z coordinate in rotated voxel space + float x_r = (r * spacing.z) - origin.z; // converted to image space + + for (int dq = -1; dq <= 1; dq++) + { + int q = (int)roundf(q_r) + dq; + float x_q = (q * spacing.y) - origin.y; + + for (int dp = -1; dp <= 1; dp++) + { + int p = (int)roundf(p_r) + dp; + float x_p = (p * spacing.x) - origin.x; + + // get distance between nearest neighbour and central voxel + // both need to be in image space + float delta_x = X_rot_x - x_p; + float delta_y = X_rot_y - x_q; + float delta_z = X_rot_z - x_r; + + // caulculate gaussian kernel + float g = expf(-1. * ((delta_x * delta_x) + (delta_y * delta_y) + (delta_z * delta_z)) / (2 * sigma * sigma)) G + += g; + }; + }; }; - }; // loop again but this time actually fetch the image values float accumulation = 0; - for (int dr = -1; dr <= 1; dr++) { - float r = dr + (int)roundf(r_r); - float x_r = origin.z + (r * spacing.z); - - for (int dq = -1; dq <= 1; dq++) { - float q = dq + (int)roundf(q_r); - float x_q = origin.y + (q * spacing.y); - - for (int dp = -1; dp <= 1; dp++) { - float p = dp + (int)roundf(p_r); - float x_p = origin.x + (p * spacing.x); - - float delta_x = Xr_x - x_p; - float delta_y = Xr_y - x_q; - float delta_z = Xr_z - x_r; - - // get input image voxel index that we are at - int i_idx = p + (q * dim.x) + (r * dim.x * dim.y); - - // Pushing the weighted counts to NN - float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma)); - atomicAdd(&out_im[i_idx], in_im[idx]*g/G); - - }; + for (int dr = -1; dr <= 1; dr++) + { + float r = dr + (int)roundf(r_r); + float x_r = origin.z + (r * spacing.z); + + for (int dq = -1; dq <= 1; dq++) + { + float q = dq + (int)roundf(q_r); + float x_q = origin.y + (q * spacing.y); + + for (int dp = -1; dp <= 1; dp++) + { + float p = dp + (int)roundf(p_r); + float x_p = origin.x + (p * spacing.x); + + float delta_x = Xr_x - x_p; + float delta_y = Xr_y - x_q; + float delta_z = Xr_z - x_r; + + // get input image voxel index that we are at + int i_idx = p + (q * dim.x) + (r * dim.x * dim.y); + + // Pushing the weighted counts to NN + float g = expf(-1. * ((delta_x * delta_x) + (delta_y * delta_y) + (delta_z * delta_z)) / (2 * sigma * sigma)); + atomicAdd(&out_im[i_idx], in_im[idx] * g / G); + }; + }; }; - }; } From aec9adf5a294ca8dce0ecc9a75f1a5f1a38347b5 Mon Sep 17 00:00:00 2001 From: danieldeidda Date: Fri, 6 Mar 2026 16:34:56 +0000 Subject: [PATCH 8/9] joining all contributions; created an actual_forward_project() that calls the kernels; we need to double check that the kernel can write into the sino but probably we need to dodevice_to_host(stir_sino,cuda_array_created_by_forward) --- .../BackProjectorByBinSPECTGPU.h | 7 + .../ForwardProjectorByBinSPECTGPU.h | 13 +- .../BackProjectorByBinSPECTGPU.cxx | 6 +- .../ForwardProjectorByBinSPECTGPU.cxx | 92 +++++++++- .../SPECTGPU_projector/SPECTGPUHelper.cxx | 159 +++++------------- .../SPECTGPU_projector/SPECTGPUProjection.cu | 59 +++---- .../SPECTGPURotateAndGaussianInterpolate.cu | 98 ++++++++++- 7 files changed, 272 insertions(+), 162 deletions(-) diff --git a/src/include/stir/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.h b/src/include/stir/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.h index 023cc75cab..085cb99418 100644 --- a/src/include/stir/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.h +++ b/src/include/stir/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.h @@ -82,6 +82,13 @@ class BackProjectorByBinSPECTGPU : public RegisteredParsingObject& stir_image, + const RelatedViewgrams&, + const int min_axial_pos_num, + const int max_axial_pos_num, + const int min_tangential_pos_num, + const int max_tangential_pos_num); + private: shared_ptr _symmetries_sptr; SPECTGPUHelper _helper; diff --git a/src/include/stir/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.h b/src/include/stir/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.h index 93868fe953..ccd01ebd50 100644 --- a/src/include/stir/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.h +++ b/src/include/stir/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.h @@ -94,8 +94,19 @@ class ForwardProjectorByBinSPECTGPU : public RegisteredParsingObject _symmetries_sptr; shared_ptr _projected_data_sptr; SPECTGPUHelper _helper; int _cuda_device; diff --git a/src/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.cxx b/src/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.cxx index e02edc95f3..29c55d7e2d 100644 --- a/src/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.cxx +++ b/src/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.cxx @@ -120,11 +120,11 @@ BackProjectorByBinSPECTGPU::start_accumulating_in_new_target() } void -BackProjectorByBinSPECTGPU::actual_back_project( +BackProjectorByBinSPECTGPU::actual_back_project(DiscretisedDensity<3, float> &stir_image, const RelatedViewgrams& related_viewgrams, const int, const int, const int, const int) { - for (stir::RelatedViewgrams::const_iterator iter = related_viewgrams.begin(); iter != related_viewgrams.end(); ++iter) - _helper.convert_viewgram_stir_to_SPECTGPU(_np_sino_w_gaps, *iter); + +// call the kernels for backward } END_NAMESPACE_STIR diff --git a/src/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.cxx b/src/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.cxx index 20431dacf9..82da38a51a 100644 --- a/src/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.cxx +++ b/src/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.cxx @@ -61,7 +61,35 @@ ForwardProjectorByBinSPECTGPU::set_up(const shared_ptr& proj { ForwardProjectorByBin::set_up(proj_data_info_sptr, density_info_sptr); check(*proj_data_info_sptr, *_density_sptr); - _symmetries_sptr.reset(new TrivialDataSymmetriesForBins(proj_data_info_sptr)); + + auto& target_cast = dynamic_cast&>(*target_sptr); + auto sizes = target_cast.get_lengths(); + + this->z_dim = sizes[1]; + this->y_dim = sizes[2]; + this->x_dim = sizes[3]; + + // Set the thread block and grid dimensions using std::tuple + this->block_dim.x = 8; + this->block_dim.y = 8; + this->block_dim.z = 8; + + this->grid_dim.x = (this->x_dim + this->block_dim.x - 1) / this->block_dim.x; + this->grid_dim.y = (this->y_dim + this->block_dim.y - 1) / this->block_dim.y; + this->grid_dim.z = (this->z_dim + this->block_dim.z - 1) / this->block_dim.z; + + // Check if z_dim is 1 or only 2D is true and return an error if it is + if (this->z_dim == 1 || this->only_2D) + { + error(" requires a 3D image and only works for a 3x3x3 neighbourhood"); + return Succeeded::no; + } + +// { +// if (this->d_kappa_data) +// cudaFree(this->d_kappa_data); +// auto kappa_ptr = this->get_kappa_sptr(); +// const bool do_kappa = !is_null_ptr(kappa_ptr); // Initialise projected_data_sptr from this->_proj_data_info_sptr @@ -78,20 +106,76 @@ ForwardProjectorByBinSPECTGPU::set_up(const shared_ptr& proj void ForwardProjectorByBinSPECTGPU::actual_forward_project( - RelatedViewgrams&, const DiscretisedDensity<3, float>&, const int, const int, const int, const int) + RelatedViewgrams& stir_sino, + const DiscretisedDensity<3, float>& stir_image, + const int min_ax, + const int max_ax, + const int min_tg, + const int max_tg) { - throw std::runtime_error("Need to use set_input() if wanting to use ForwardProjectorByBinSPECTGPU."); + //for all views in relateViewgram call the kernels + dim3 cuda_block_dim(this->block_dim.x, this->block_dim.y, this->block_dim.z); + dim3 cuda_grid_dim(this->grid_dim.x, this->grid_dim.y, this->grid_dim.z); + viewgrams = _projected_data_sptr->get_related_viewgrams(viewgrams.get_basic_view_segment_num(), _symmetries_sptr); + + for (auto view=0; view min_ind, max_ind; + + stir_image.get_regular_range(min_ind, max_ind); + + const int min_z = min_ind[1]; + const int max_z = max_ind[1]; + + const int min_y = min_ind[2]; + const int max_y = max_ind[2]; + + const int min_x = min_ind[3]; + const int max_x = max_ind[3]; + + int3 dim(max_x - min_x + 1, max_y - min_y + 1, max_z - min_z + 1 ); + + float angle_rad; //todo + float3 spacing(,,); + float3 origin(,,); + + rotateKernel_pull<<>>(dev_image, + out_im, + dim); + + forwardKernel<<>>(out_im, + viewgrams[view], + dim); + + } + // cudaMalloc(&this->cuda_image, stir_image_sptr->size_all() * sizeof(elemT)); + // array_to_device(this->cuda_image, *stir_image_sptr); + } + + } void ForwardProjectorByBinSPECTGPU::actual_forward_project( RelatedViewgrams& viewgrams, const int, const int, const int, const int) { - // if (min_axial_pos_num != _proj_data_info_sptr->get_min_axial_pos_num() || + if (min_axial_pos_num != _proj_data_info_sptr->get_min_axial_pos_num() || // ... ) // error(); +//for all views in relateViewgram call the kernels viewgrams = _projected_data_sptr->get_related_viewgrams(viewgrams.get_basic_view_segment_num(), _symmetries_sptr); + for (auto view=0; viewcuda_image, stir_image_sptr->size_all() * sizeof(elemT)); + } + // cudaMalloc(&this->cuda_image, stir_image_sptr->size_all() * sizeof(elemT)); +// array_to_device(this->cuda_image, *stir_image_sptr); } void diff --git a/src/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.cxx b/src/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.cxx index 2ca9788968..1325576b09 100644 --- a/src/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.cxx +++ b/src/recon_buildblock/SPECTGPU_projector/SPECTGPUHelper.cxx @@ -31,91 +31,57 @@ #include "stir/IO/stir_ecat_common.h" #include "stir/error.h" #include "stir/format.h" +#include "stir/cuda_utilities.h" // Non-STIR includes #include #include #include "driver_types.h" // SPECTGPU includes -#include "def.h" -#include "auxmath.h" -#include "prjb.h" -#include "prjf.h" -#include "recon.h" -#include "lmproc.h" -#include "scanner_0.h" -#include "rnd.h" -#include "norm.h" START_NAMESPACE_STIR SPECTGPUHelper::~SPECTGPUHelper() {} -// static void -// delete_axialLUT(axialLUT* axlut_ptr) -// { -// if (!axlut_ptr) -// return; -// delete[] axlut_ptr->li2rno; -// delete[] axlut_ptr->li2sn; -// delete[] axlut_ptr->li2nos; -// delete[] axlut_ptr->sn1_rno; -// delete[] axlut_ptr->sn1_sn11; -// delete[] axlut_ptr->sn1_ssrb; -// delete[] axlut_ptr->sn1_sn11no; -// } - -static shared_ptr -get_cnst(const Scanner& scanner, const bool cuda_verbose, const char cuda_device) -{ - shared_ptr cnt_sptr = MAKE_SHARED(); - cnt_sptr->DEVID = cuda_device; // device (GPU) ID. allows choosing the device on which to perform calculations - cnt_sptr->VERBOSE = cuda_verbose; +//static shared_ptr +//get_cnst(const Scanner& scanner, const bool cuda_verbose, const char cuda_device) +//{ +// shared_ptr cnt_sptr = MAKE_SHARED(); + +// cnt_sptr->DEVID = cuda_device; // device (GPU) ID. allows choosing the device on which to perform calculations +// cnt_sptr->VERBOSE = cuda_verbose; - cnt_sptr->A = NSANGLES; // sino angles - cnt_sptr->W = NSBINS; // sino bins for any angular index - cnt_sptr->aw = AW; // sino bins (active only) +// cnt_sptr->A = NSANGLES; // sino angles +// cnt_sptr->W = NSBINS; // sino bins for any angular index +// cnt_sptr->aw = AW; // sino bins (active only) - cnt_sptr->NCRS = nCRS; // number of crystals - cnt_sptr->NRNG = NRINGS; // number of axial positions - cnt_sptr->D = -1; // number of linear indexes along Michelogram diagonals /*unknown*/ - cnt_sptr->Bt = -1; // number of buckets transaxially /*unknown*/ +// cnt_sptr->NCRS = nCRS; // number of crystals +// cnt_sptr->NRNG = NRINGS; // number of axial positions +// cnt_sptr->D = -1; // number of linear indexes along Michelogram diagonals /*unknown*/ +// cnt_sptr->Bt = -1; // number of buckets transaxially /*unknown*/ - cnt_sptr->B = NBUCKTS; // number of buckets (total) - cnt_sptr->Cbt = 32552; // number of crystals in bucket transaxially /*unknown*/ - cnt_sptr->Cba = 3; // number of crystals in bucket axially /*unknown*/ +// cnt_sptr->B = NBUCKTS; // number of buckets (total) +// cnt_sptr->Cbt = 32552; // number of crystals in bucket transaxially /*unknown*/ +// cnt_sptr->Cba = 3; // number of crystals in bucket axially /*unknown*/ - cnt_sptr->NSN1 = NSINOS; // number of sinos - cnt_sptr->NSN64 = NRINGS * NRINGS; // with no MRD limit - cnt_sptr->NSEG0 = SEG0; +// cnt_sptr->NSN1 = NSINOS; // number of sinos +// cnt_sptr->NSN64 = NRINGS * NRINGS; // with no MRD limit +// cnt_sptr->NSEG0 = SEG0; - cnt_sptr->RNG_STRT = 0; - cnt_sptr->RNG_END = NRINGS; +// cnt_sptr->RNG_STRT = 0; +// cnt_sptr->RNG_END = NRINGS; - cnt_sptr->ALPHA = aLPHA; // angle subtended by a crystal - float R = 32.8f; // ring radius - cnt_sptr->RE = R + 0.67f; // effective ring radius accounting for the depth of interaction - cnt_sptr->AXR = SZ_RING; // axial crystal dim +// cnt_sptr->ALPHA = aLPHA; // angle subtended by a crystal +// float R = 32.8f; // ring radius +// cnt_sptr->RE = R + 0.67f; // effective ring radius accounting for the depth of interaction +// cnt_sptr->AXR = SZ_RING; // axial crystal dim - float CLGHT = 29979245800.f; // speed of light [cm/s] - return cnt_sptr; -} +// float CLGHT = 29979245800.f; // speed of light [cm/s] +// return cnt_sptr; +//} -static inline unsigned -to_1d_idx(const unsigned nrow, const unsigned ncol, const unsigned row, const unsigned col) -{ - return col + ncol * row; -} -template -dataType* -create_heap_array(const unsigned numel, const dataType val = dataType(0)) -{ - dataType* array = new dataType[numel]; - std::fill(array, array + numel, val); - return array; -} void SPECTGPUHelper::set_up() @@ -124,14 +90,14 @@ SPECTGPUHelper::set_up() throw std::runtime_error("SPECTGPUHelper::set_up() " "emission or transmission mode (att) not set."); - // Get consts - _cnt_sptr = get_cnst(_scanner_type, _verbose, _devid); +// // Get consts +// _cnt_sptr = get_cnst(_scanner_type, _verbose, _devid); - // isub - _isub = std::vector(unsigned(AW)); - for (unsigned i = 0; i < unsigned(AW); i++) - _isub[i] = int(i); +// // isub +// _isub = std::vector(unsigned(AW)); +// for (unsigned i = 0; i < unsigned(AW); i++) +// _isub[i] = int(i); _already_set_up = true; } @@ -196,53 +162,6 @@ SPECTGPUHelper::convert_SPECTGPU_proj_3d_to_1d_idx(const unsigned ang, const uns } void -SPECTGPUHelper::permute(std::vector& output_array, - const std::vector& orig_array, - const unsigned output_dims[3], - const unsigned permute_order[3]) const -{ -#ifndef NDEBUG - // Check that in the permute order, each number is between 0 and 2 (can't be <0 because it's unsigned) - for (unsigned i = 0; i < 3; ++i) - if (permute_order[i] > 2) - throw std::runtime_error("Permute order values should be between 0 and 2."); - // Check that each number is unique - for (unsigned i = 0; i < 3; ++i) - for (unsigned j = i + 1; j < 3; ++j) - if (permute_order[i] == permute_order[j]) - throw std::runtime_error("Permute order values should be unique."); - // Check that size of output_dims==arr.size() - assert(orig_array.size() == output_dims[0] * output_dims[1] * output_dims[2]); - // Check that output array is same size as input array - assert(orig_array.size() == output_array.size()); -#endif - - // Calculate old dimensions - unsigned old_dims[3]; - for (unsigned i = 0; i < 3; ++i) - old_dims[permute_order[i]] = output_dims[i]; - - // Loop over all elements - for (unsigned old_1d_idx = 0; old_1d_idx < orig_array.size(); ++old_1d_idx) - { - - // From the 1d index, generate the old 3d index - unsigned old_3d_idx[3] - = { old_1d_idx / (old_dims[2] * old_dims[1]), (old_1d_idx / old_dims[2]) % old_dims[1], old_1d_idx % old_dims[2] }; - - // Get the corresponding new 3d index - unsigned new_3d_idx[3]; - for (unsigned i = 0; i < 3; ++i) - new_3d_idx[i] = old_3d_idx[permute_order[i]]; - - // Get the new 1d index from the new 3d index - const unsigned new_1d_idx - = new_3d_idx[0] * output_dims[2] * output_dims[1] + new_3d_idx[1] * output_dims[2] + new_3d_idx[2]; - - // Fill the data - output_array[new_1d_idx] = orig_array[old_1d_idx]; - } -} void SPECTGPUHelper::back_project(std::vector& image, const std::vector& sino_no_gaps) const @@ -262,11 +181,13 @@ SPECTGPUHelper::back_project(std::vector& image, const std::vector } void -SPECTGPUHelper::forward_project(std::vector& sino, const std::vector& image) const +SPECTGPUHelper::forward_project(Array<3, elemT>& sino, const Array<3, elemT>& image) const { check_set_up(); assert(!sino.empty()); - +// prjdatainmemory + cudaMalloc(&this->cuda_image, stir_image_sptr->size_all() * sizeof(elemT)); + array_to_device(this->cuda_image, *stir_image_sptr); // Permute the data (as this is done on the SPECTGPU python side before forward projection // unsigned output_dims[3] = { 320, 320, 127 }; // unsigned permute_order[3] = { 1, 2, 0 }; diff --git a/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu b/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu index 6aa040140b..b6fbbdc779 100644 --- a/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu +++ b/src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu @@ -26,44 +26,37 @@ #include #include - //the following is a pull operation __global__ void forwardKernel(const float* __restrict__ in_im, float* __restrict__ out_sino, - int3 dim, - float3 spacing, - float3 origin, - float angle_rad) - { - //1. define position in the detector, x horizontal, z vertical, y is the depth. - - int det_x = blockIdx.x * blockDim.x + threadIdx.x; - int det_y = blockIdx.y * blockDim.y + threadIdx.y; - int det_z = blockIdx.z * blockDim.z + threadIdx.z; - - if (det_x >= dim.x || det_y >= dim.y || det_z >= dim.z) { - return; - } - int sino_idx = det_z * dim.x + det_x; + int3 dim) +{ + //1. define position in the detector, x horizontal, z vertical, y is the depth. + + int det_x = blockIdx.x * blockDim.x + threadIdx.x; + int det_y = blockIdx.y * blockDim.y + threadIdx.y; + int det_z = blockIdx.z * blockDim.z + threadIdx.z; + + if (det_x >= dim.x || det_y >= dim.y || det_z >= dim.z) { + return; + } + int sino_idx = det_z * dim.x + det_x; + + // Sum voxel values along the y-axis (depth) for this detector pixel + int voxel_idx = det_y * dim.x * dim.z + det_z * dim.x + det_x; + out_sino[sino_idx] += in_im[voxel_idx]; - // Sum voxel values along the y-axis (depth) for this detector pixel - int voxel_idx = det_y * dim.x * dim.z + det_z * dim.x + det_x; - out_sino[sino_idx] += in_im[voxel_idx]; - - } +} // the following is the adjoint operation (push) __global__ void backwardKernel(const float* __restrict__ in_sino, float* __restrict__ out_im, - int3 dim, - float3 spacing, - float3 origin, - float angle_rad) - { - int idz = blockIdx.zblockDim.z + threadIdx.z; - int idy = blockIdx.yblockDim.y + threadIdx.y; - int idx = blockIdx.xblockDim.x + threadIdx.x; - if ( idz < dim.z && idy < dim.y && idx < dim.x) - { - out_im[(idz dim.y + idy) * dim.x + idx] += sinogram[idz * dim.y + idy] * spacing.x; - } + int3 dim) +{ + int idz = blockIdx.zblockDim.z + threadIdx.z; + int idy = blockIdx.yblockDim.y + threadIdx.y; + int idx = blockIdx.xblockDim.x + threadIdx.x; + if ( idz < dim.z && idy < dim.y && idx < dim.x) + { + out_im[(idz dim.y + idy) * dim.x + idx] += sinogram[idz * dim.y + idy] * spacing.x; + } } diff --git a/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu b/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu index 1e19382d67..7d425ea368 100644 --- a/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu +++ b/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu @@ -136,6 +136,7 @@ rotateKernel_pull( }; // the following is the adjoint operation (push) +<<<<<<< Updated upstream __global__ void rotateKernel_push( const float* __restrict__ in_im, float* __restrict__ out_im, int3 dim, float3 spacing, float3 origin, float angle_rad) @@ -202,13 +203,81 @@ rotateKernel_push( // caulculate gaussian kernel float g = expf(-1. * ((delta_x * delta_x) + (delta_y * delta_y) + (delta_z * delta_z)) / (2 * sigma * sigma)) G += g; +======= +__global__ void rotateKernel_push(const float* __restrict__ in_im, + float* __restrict__ out_im, + int3 dim, + float3 spacing, + float3 origin, + float angle_rad) +{ + // parallelise the operation across all image voxels + int i = threadIdx.x + blockDim.x * blockIdx.x; + int j = threadIdx.y + blockDim.y * blockIdx.y; + int k = threadIdx.z + blockDim.z * blockIdx.z; + + // check we are not outside the image + if (i >= dim.x || j >= dim.y || k >= dim.z) return; + + // get the 1-dimensional index + int idx = i + (j * dim.x) + (k * dim.x * dim.y); + + // move to centre of voxel and convert from voxel to image space + float Xs_rel_x = (i + 0.5f) * spacing.x - origin.x; + float Xs_rel_y = (j + 0.5f) * spacing.y - origin.y; + float Xs_rel_z = (k + 0.5f) * spacing.z - origin.z; + + // perform the rotation by 'angle_rad' + float X_rot_x = (Xs_rel_x * cosf(angle_rad)) - (Xs_rel_y * sinf(angle_rad)); + float X_rot_y = (Xs_rel_x * sinf(angle_rad)) + (Xs_rel_y * cosf(angle_rad)); + float X_rot_z = Xs_rel_z; + + // convert back to voxel space origin (corner of image) but keep image space dimensions + float Xr_x = X_rot_x + origin.x; + float Xr_y = X_rot_y + origin.y; + float Xr_z = X_rot_z + origin.z; + + // now we go fully back to voxel space (but now it's been rotated) + float p_r = Xr_x / spacing.x; + float q_r = Xr_y / spacing.y; + float r_r = Xr_z / spacing.z; + + // first loop is for finding the value of a gaussian kernel at + // each nearest neighbour position and summing them all. This will allow normalisation + // of each NN contribution in the next loop. + float G = 0; // accumulator variable + float sigma = 1; // gaussian kernel sigma + + for (int dr = -1; dr <= 1; dr++) { + int r = (int)roundf(r_r) + dr; // nearest neighbour z coordinate in rotated voxel space + float x_r = (r * spacing.z) - origin.z; // converted to image space + + for (int dq = -1; dq <= 1; dq++) { + int q = (int)roundf(q_r) + dq; + float x_q = (q * spacing.y) - origin.y; + + for (int dp = -1; dp <= 1; dp++) { + int p = (int)roundf(p_r) + dp; + float x_p = (p * spacing.x) - origin.x; + + // get distance between nearest neighbour and central voxel + // both need to be in image space + float delta_x = X_rot_x - x_p; + float delta_y = X_rot_y - x_q; + float delta_z = X_rot_z - x_r; + + // caulculate gaussian kernel + float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma)) + G += g; +>>>>>>> Stashed changes }; }; }; - // loop again but this time actually fetch the image values - float accumulation = 0; + // loop again but this time actually fetch the image values + float accumulation = 0; +<<<<<<< Updated upstream for (int dr = -1; dr <= 1; dr++) { float r = dr + (int)roundf(r_r); @@ -234,6 +303,31 @@ rotateKernel_push( // Pushing the weighted counts to NN float g = expf(-1. * ((delta_x * delta_x) + (delta_y * delta_y) + (delta_z * delta_z)) / (2 * sigma * sigma)); atomicAdd(&out_im[i_idx], in_im[idx] * g / G); +======= + for (int dr = -1; dr <= 1; dr++) { + float r = dr + (int)roundf(r_r); + float x_r = origin.z + (r * spacing.z); + + for (int dq = -1; dq <= 1; dq++) { + float q = dq + (int)roundf(q_r); + float x_q = origin.y + (q * spacing.y); + + for (int dp = -1; dp <= 1; dp++) { + float p = dp + (int)roundf(p_r); + float x_p = origin.x + (p * spacing.x); + + float delta_x = Xr_x - x_p; + float delta_y = Xr_y - x_q; + float delta_z = Xr_z - x_r; + + // get input image voxel index that we are at + int i_idx = p + (q * dim.x) + (r * dim.x * dim.y); + + // Pushing the weighted counts to NN + float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma)); + atomicAdd(&out_im[i_idx], in_im[idx]*g/G); + +>>>>>>> Stashed changes }; }; }; From edc60a858d19b01a7b90a61f21a5e7ea7e2d0ce7 Mon Sep 17 00:00:00 2001 From: "david.roddy" Date: Fri, 6 Mar 2026 16:49:18 +0000 Subject: [PATCH 9/9] fixup! joining all contributions; created an actual_forward_project() that calls the kernels; we need to double check that the kernel can write into the sino but probably we need to dodevice_to_host(stir_sino,cuda_array_created_by_forward) --- .../SPECTGPURotateAndGaussianInterpolate.cu | 98 +------------------ 1 file changed, 2 insertions(+), 96 deletions(-) diff --git a/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu b/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu index 7d425ea368..1e19382d67 100644 --- a/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu +++ b/src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu @@ -136,7 +136,6 @@ rotateKernel_pull( }; // the following is the adjoint operation (push) -<<<<<<< Updated upstream __global__ void rotateKernel_push( const float* __restrict__ in_im, float* __restrict__ out_im, int3 dim, float3 spacing, float3 origin, float angle_rad) @@ -203,81 +202,13 @@ rotateKernel_push( // caulculate gaussian kernel float g = expf(-1. * ((delta_x * delta_x) + (delta_y * delta_y) + (delta_z * delta_z)) / (2 * sigma * sigma)) G += g; -======= -__global__ void rotateKernel_push(const float* __restrict__ in_im, - float* __restrict__ out_im, - int3 dim, - float3 spacing, - float3 origin, - float angle_rad) -{ - // parallelise the operation across all image voxels - int i = threadIdx.x + blockDim.x * blockIdx.x; - int j = threadIdx.y + blockDim.y * blockIdx.y; - int k = threadIdx.z + blockDim.z * blockIdx.z; - - // check we are not outside the image - if (i >= dim.x || j >= dim.y || k >= dim.z) return; - - // get the 1-dimensional index - int idx = i + (j * dim.x) + (k * dim.x * dim.y); - - // move to centre of voxel and convert from voxel to image space - float Xs_rel_x = (i + 0.5f) * spacing.x - origin.x; - float Xs_rel_y = (j + 0.5f) * spacing.y - origin.y; - float Xs_rel_z = (k + 0.5f) * spacing.z - origin.z; - - // perform the rotation by 'angle_rad' - float X_rot_x = (Xs_rel_x * cosf(angle_rad)) - (Xs_rel_y * sinf(angle_rad)); - float X_rot_y = (Xs_rel_x * sinf(angle_rad)) + (Xs_rel_y * cosf(angle_rad)); - float X_rot_z = Xs_rel_z; - - // convert back to voxel space origin (corner of image) but keep image space dimensions - float Xr_x = X_rot_x + origin.x; - float Xr_y = X_rot_y + origin.y; - float Xr_z = X_rot_z + origin.z; - - // now we go fully back to voxel space (but now it's been rotated) - float p_r = Xr_x / spacing.x; - float q_r = Xr_y / spacing.y; - float r_r = Xr_z / spacing.z; - - // first loop is for finding the value of a gaussian kernel at - // each nearest neighbour position and summing them all. This will allow normalisation - // of each NN contribution in the next loop. - float G = 0; // accumulator variable - float sigma = 1; // gaussian kernel sigma - - for (int dr = -1; dr <= 1; dr++) { - int r = (int)roundf(r_r) + dr; // nearest neighbour z coordinate in rotated voxel space - float x_r = (r * spacing.z) - origin.z; // converted to image space - - for (int dq = -1; dq <= 1; dq++) { - int q = (int)roundf(q_r) + dq; - float x_q = (q * spacing.y) - origin.y; - - for (int dp = -1; dp <= 1; dp++) { - int p = (int)roundf(p_r) + dp; - float x_p = (p * spacing.x) - origin.x; - - // get distance between nearest neighbour and central voxel - // both need to be in image space - float delta_x = X_rot_x - x_p; - float delta_y = X_rot_y - x_q; - float delta_z = X_rot_z - x_r; - - // caulculate gaussian kernel - float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma)) - G += g; ->>>>>>> Stashed changes }; }; }; - // loop again but this time actually fetch the image values - float accumulation = 0; + // loop again but this time actually fetch the image values + float accumulation = 0; -<<<<<<< Updated upstream for (int dr = -1; dr <= 1; dr++) { float r = dr + (int)roundf(r_r); @@ -303,31 +234,6 @@ __global__ void rotateKernel_push(const float* __restrict__ in_im, // Pushing the weighted counts to NN float g = expf(-1. * ((delta_x * delta_x) + (delta_y * delta_y) + (delta_z * delta_z)) / (2 * sigma * sigma)); atomicAdd(&out_im[i_idx], in_im[idx] * g / G); -======= - for (int dr = -1; dr <= 1; dr++) { - float r = dr + (int)roundf(r_r); - float x_r = origin.z + (r * spacing.z); - - for (int dq = -1; dq <= 1; dq++) { - float q = dq + (int)roundf(q_r); - float x_q = origin.y + (q * spacing.y); - - for (int dp = -1; dp <= 1; dp++) { - float p = dp + (int)roundf(p_r); - float x_p = origin.x + (p * spacing.x); - - float delta_x = Xr_x - x_p; - float delta_y = Xr_y - x_q; - float delta_z = Xr_z - x_r; - - // get input image voxel index that we are at - int i_idx = p + (q * dim.x) + (r * dim.x * dim.y); - - // Pushing the weighted counts to NN - float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma)); - atomicAdd(&out_im[i_idx], in_im[idx]*g/G); - ->>>>>>> Stashed changes }; }; };