From 94efb22e2a2e1d6974e282a6c0d0ca6012eb667f Mon Sep 17 00:00:00 2001 From: nmdicom-recon Date: Thu, 2 Mar 2023 13:20:59 +0000 Subject: [PATCH 1/4] making useful global variable class members --- .../recon_buildblock/ProjMatrixByBinSPECTUB.h | 13 ++++--- .../stir/recon_buildblock/SPECTUB_Weight3d.h | 29 ++++++++------- .../ProjMatrixByBinSPECTUB.cxx | 29 +++++++++------ src/recon_buildblock/SPECTUB_Weight3d.cxx | 35 ++++++++++--------- 4 files changed, 63 insertions(+), 43 deletions(-) diff --git a/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h b/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h index eae1ecd495..51022e6ada 100644 --- a/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h +++ b/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h @@ -75,7 +75,7 @@ class Bin; End Projection Matrix By Bin SPECT UB Parameters:= \endverbatim */ - +using namespace SPECTUB; class ProjMatrixByBinSPECTUB : public RegisteredParsingObject< ProjMatrixByBinSPECTUB, @@ -97,7 +97,7 @@ class ProjMatrixByBinSPECTUB : virtual void set_up( const shared_ptr& proj_data_info_ptr, const shared_ptr >& density_info_ptr // TODO should be Info only - ); +); bool get_keep_all_views_in_cache() const; //! Enable keeping the matrix in memory @@ -146,7 +146,6 @@ class ProjMatrixByBinSPECTUB : void set_resolution_model(const float collimator_sigma_0_in_mm, const float collimator_slope_in_mm, const bool full_3D = true); - private: // parameters that will be parsed @@ -172,6 +171,9 @@ class ProjMatrixByBinSPECTUB : bool already_setup; + static wm_da_type wm; + static wmh_type wmh; // this could be an arry of wmh_type for each index + float * Rrad; virtual void calculate_proj_matrix_elems_for_one_bin( @@ -206,7 +208,10 @@ class ProjMatrixByBinSPECTUB : int maxszb; - void compute_one_subset(const int kOS) const; + void compute_one_subset(const int kOS, + wm_da_type& wm, + wmh_type& wmh, + float * Rrad) const; void delete_UB_SPECT_arrays(); mutable std::vector subset_already_processed; }; diff --git a/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h b/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h index 5b5f28d88d..bdc0848a49 100644 --- a/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h +++ b/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h @@ -51,20 +51,25 @@ void wm_calculation( const int kOS, const bool *msk_2d, const int maxszb, const discrf_type *const gaussdens, - const int *const NITEMS + const int *const NITEMS, + wm_da_type& wm, + wmh_type& wmh, + float * Rrad ); void wm_size_estimation (int kOS, - const angle_type * const ang, - voxel_type vox, - bin_type bin, - const volume_type& vol, - const proj_type& prj, - const bool * const msk_3d, - const bool *const msk_2d, - const int maxszb, - const discrf_type * const gaussdens, - int *NITEMS); + const angle_type * const ang, + voxel_type vox, + bin_type bin, + const volume_type& vol, + const proj_type& prj, + const bool * const msk_3d, + const bool *const msk_2d, + const int maxszb, + const discrf_type * const gaussdens, + int *NITEMS, + wmh_type &wmh, + float *Rrad); //... geometric component ............................................ @@ -96,7 +101,7 @@ void calc_psf_bin ( float center_psf, float binszcm, discrf_type const * const v void calc_att_path ( const bin_type& bin, const voxel_type& vox, const volume_type& vol, attpth_type *attpth ); -float calc_att ( const attpth_type *const attpth, const float *const attmap, int islc ); +float calc_att (const attpth_type *const attpth, const float *const attmap, int islc , wmh_type &wmh); int comp_dist ( float dx, float dy, float dz, float dlast ); diff --git a/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx b/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx index c6ae155aa2..0164629f17 100644 --- a/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx +++ b/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx @@ -2,6 +2,7 @@ Copyright (C) 2013, Institute for Bioengineering of Catalonia Copyright (C) Biomedical Image Group (GIB), Universitat de Barcelona, Barcelona, Spain. Copyright (C) 2013-2014, 2019, 2020 University College London + Copyright (C) 2023 National Physical Laboratory This file is part of STIR. SPDX-License-Identifier: Apache-2.0 @@ -17,6 +18,7 @@ \author Berta Marti Fuster \author Carles Falcon \author Kris Thielemans + \author Daniel Deidda */ //#include "stir/ProjDataInterfile.h" @@ -60,14 +62,14 @@ #include "stir/recon_buildblock/SPECTUB_Weight3d.h" /* UB-SPECT global variables */ -namespace SPECTUB { - wm_da_type wm; - wmh_type wmh; - float * Rrad; -} +//namespace SPECTUB { +// wm_da_type wm; +// wmh_type wmh; +// float * Rrad; +//} START_NAMESPACE_STIR - +using namespace SPECTUB; const char * const ProjMatrixByBinSPECTUB::registered_name = @@ -658,7 +660,8 @@ set_up( //... size estimations ........................................................ - wm_size_estimation ( kOS, ang, vox, bin, vol, prj, msk_3d, msk_2d, maxszb, &gaussdens, NITEMS[kOS] ); + wm_size_estimation ( kOS, ang, vox, bin, vol, prj, msk_3d, msk_2d, maxszb, &gaussdens, NITEMS[kOS], + wmh,Rrad ); //cout << "\nwm_SPECT. Size estimation done. time (s): " << double( clock()-ini )/CLOCKS_PER_SEC < Date: Fri, 10 Mar 2023 11:00:15 +0000 Subject: [PATCH 2/4] making variable non-static; removing unnecessary namespace; --- .../recon_buildblock/ProjMatrixByBinSPECTUB.h | 11 ++-- .../stir/recon_buildblock/SPECTUB_Weight3d.h | 64 +++++++++---------- .../ProjMatrixByBinSPECTUB.cxx | 56 ++++++++-------- src/recon_buildblock/SPECTUB_Weight3d.cxx | 49 ++++++++------ 4 files changed, 89 insertions(+), 91 deletions(-) diff --git a/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h b/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h index 51022e6ada..4127343d01 100644 --- a/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h +++ b/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h @@ -75,7 +75,7 @@ class Bin; End Projection Matrix By Bin SPECT UB Parameters:= \endverbatim */ -using namespace SPECTUB; +//using namespace SPECTUB; class ProjMatrixByBinSPECTUB : public RegisteredParsingObject< ProjMatrixByBinSPECTUB, @@ -90,6 +90,9 @@ class ProjMatrixByBinSPECTUB : //! Default constructor (calls set_defaults()) ProjMatrixByBinSPECTUB(); + // disable copy-constructor as currently unsafe to copy due to bare pointers + ProjMatrixByBinSPECTUB(const ProjMatrixByBinSPECTUB&) = delete; + //! Destructor (deallocates UB SPECT memory) ~ProjMatrixByBinSPECTUB(); @@ -171,8 +174,8 @@ class ProjMatrixByBinSPECTUB : bool already_setup; - static wm_da_type wm; - static wmh_type wmh; // this could be an arry of wmh_type for each index + mutable SPECTUB::wm_da_type wm; + mutable SPECTUB::wmh_type wmh; // this could be an arry of wmh_type for each index float * Rrad; virtual void @@ -209,8 +212,6 @@ class ProjMatrixByBinSPECTUB : void compute_one_subset(const int kOS, - wm_da_type& wm, - wmh_type& wmh, float * Rrad) const; void delete_UB_SPECT_arrays(); mutable std::vector subset_already_processed; diff --git a/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h b/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h index bdc0848a49..57888f662f 100644 --- a/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h +++ b/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h @@ -23,73 +23,67 @@ Berta Marti Fuster, Carles Falcon, Charalampos Tsoumpas, Lefteris Livieratos, Pablo Aguiar, Albert Cot, Domenec Ros and Kris Thielemans, Integration of advanced 3D SPECT modeling into the open-source STIR framework, Med. Phys. 40, 092502 (2013); http://dx.doi.org/10.1118/1.4816676 - - \todo Variables wm, wmh and Rrad are currently global variables. This means that this code would be very dangerous - in a parallel setting. */ namespace SPECTUB { -/* Global variables: the matrix etc TODO - Need to be defined elsewhere - - A lot of the functions below modify these variables. -*/ - extern wm_da_type wm; //! weight (or probability) matrix - extern wmh_type wmh; //! information to construct wm - extern float * Rrad; //! radii per view +// extern wm_da_type wm; //! weight (or probability) matrix +// extern wmh_type wmh; //! information to construct wm +// extern float * Rrad; //! radii per view void wm_calculation( const int kOS, - const angle_type *const ang, - voxel_type vox, + const SPECTUB::angle_type *const ang, + SPECTUB::voxel_type vox, bin_type bin, - const volume_type& vol, + const SPECTUB::volume_type& vol, const proj_type& prj, const float *attmap, const bool *msk_3d, const bool *msk_2d, const int maxszb, - const discrf_type *const gaussdens, + const SPECTUB::discrf_type *const gaussdens, const int *const NITEMS, - wm_da_type& wm, - wmh_type& wmh, + SPECTUB::wm_da_type& wm, + SPECTUB::wmh_type& wmh, float * Rrad ); void wm_size_estimation (int kOS, - const angle_type * const ang, - voxel_type vox, + const SPECTUB::angle_type * const ang, + SPECTUB::voxel_type vox, bin_type bin, - const volume_type& vol, - const proj_type& prj, + const SPECTUB::volume_type& vol, + const SPECTUB::proj_type& prj, const bool * const msk_3d, const bool *const msk_2d, const int maxszb, - const discrf_type * const gaussdens, + const SPECTUB::discrf_type * const gaussdens, int *NITEMS, - wmh_type &wmh, + SPECTUB::wmh_type &wmh, float *Rrad); //... geometric component ............................................ -void calc_gauss ( discrf_type *gaussdens ); +void calc_gauss ( SPECTUB::discrf_type *gaussdens ); -void calc_vxprj ( angle_type *ang ); +void calc_vxprj ( SPECTUB::angle_type *ang ); -void voxel_projection ( voxel_type *vox, float * eff, float lngcmd2 ); +void voxel_projection (SPECTUB::voxel_type *vox, float * eff, float lngcmd2 , SPECTUB::wmh_type &wmh); -void fill_psf_no( psf2da_type *psf, psf1d_type *psf1d_h, const voxel_type& vox, const angle_type * const ang, float szdx ); +void fill_psf_no(SPECTUB::psf2da_type *psf, SPECTUB::psf1d_type *psf1d_h, const SPECTUB::voxel_type& vox, const angle_type * const ang, float szdx , SPECTUB::wmh_type &wmh); -void fill_psf_2d( psf2da_type *psf, psf1d_type *psf1d_h, const voxel_type& vox, discrf_type const*const gaussdens, float szdx ); +void fill_psf_2d( SPECTUB::psf2da_type *psf, SPECTUB::psf1d_type *psf1d_h, const SPECTUB::voxel_type& vox, SPECTUB::discrf_type const*const gaussdens, float szdx, SPECTUB::wmh_type &wmh); -void fill_psf_3d(psf2da_type *psf, - psf1d_type *psf1d_h, - psf1d_type *psf1d_v, - const voxel_type& vox, discrf_type const * const gaussdens, float szdx, float thdx, float thcmd2 ); +void fill_psf_3d(SPECTUB::psf2da_type *psf, + SPECTUB::psf1d_type *psf1d_h, + SPECTUB::psf1d_type *psf1d_v, + const SPECTUB::voxel_type& vox, SPECTUB::discrf_type const * const gaussdens, + float szdx, float thdx, float thcmd2, + wmh_type &wmh); -void calc_psf_bin ( float center_psf, float binszcm, discrf_type const * const vxprj, psf1d_type *psf ); +void calc_psf_bin (float center_psf, float binszcm, SPECTUB::discrf_type const * const vxprj, SPECTUB::psf1d_type *psf , SPECTUB::wmh_type &wmh); //... attenuation................................................... @@ -99,9 +93,9 @@ void calc_psf_bin ( float center_psf, float binszcm, discrf_type const * const v //void size_attpth_full( psf2da_type *psf, voxel_type vox, volume_type vol, float *att, angle_type *ang ); -void calc_att_path ( const bin_type& bin, const voxel_type& vox, const volume_type& vol, attpth_type *attpth ); +void calc_att_path ( const bin_type& bin, const SPECTUB::voxel_type& vox, const SPECTUB::volume_type& vol, SPECTUB::attpth_type *attpth ); -float calc_att (const attpth_type *const attpth, const float *const attmap, int islc , wmh_type &wmh); +float calc_att (const SPECTUB::attpth_type *const attpth, const float *const attmap, int islc , SPECTUB::wmh_type &wmh); int comp_dist ( float dx, float dy, float dz, float dlast ); diff --git a/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx b/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx index 0164629f17..79f9066a56 100644 --- a/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx +++ b/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx @@ -68,8 +68,8 @@ // float * Rrad; //} -START_NAMESPACE_STIR using namespace SPECTUB; +START_NAMESPACE_STIR const char * const ProjMatrixByBinSPECTUB::registered_name = @@ -250,8 +250,6 @@ set_up( } #endif - using namespace SPECTUB; - const VoxelsOnCartesianGrid * image_info_ptr = dynamic_cast*> (density_info_ptr.get()); @@ -689,7 +687,6 @@ delete_UB_SPECT_arrays() if (!this->already_setup) return; //... freeing matrix memory.................................... - using namespace SPECTUB; delete [] Rrad; if ( !wmh.do_psf ){ @@ -738,11 +735,8 @@ delete_UB_SPECT_arrays() void ProjMatrixByBinSPECTUB:: compute_one_subset(const int kOS, - wm_da_type& wm, - wmh_type& wmh, float * Rrad) const { - using namespace SPECTUB; CPUTimer timer; timer.start(); @@ -750,12 +744,12 @@ compute_one_subset(const int kOS, //... to fill wmh fields related to the subset .................................. - wmh.subset_ind = kOS; + this->wmh.subset_ind = kOS; for ( int i = 0 ; i < prj.NangOS ; i ++ ){ - wmh.index[ i ] = prj.order[ i + kOS * prj.NangOS ]; - wmh.Rrad [ i ] = Rrad[ wmh.index[ i ] ]; + this->wmh.index[ i ] = prj.order[ i + kOS * prj.NangOS ]; + this->wmh.Rrad [ i ] = Rrad[ this->wmh.index[ i ] ]; } //... NITEMS initialization ...................... @@ -771,26 +765,26 @@ compute_one_subset(const int kOS, int ne = 0; - for ( int i = 0 ; i < wmh.prj.NbOS ; i++ ) ne += NITEMS[kOS][ i ]; + for ( int i = 0 ; i < this->wmh.prj.NbOS ; i++ ) ne += NITEMS[kOS][ i ]; //... size information .................................................................... info(boost::format("total number of non-zero weights in this view: %1%, estimated size: %2% MB") % ne - % ( wm.do_save_STIR ? (ne + 10* prj.NbOS)/104857.6 : ne/131072), + % ( this->wm.do_save_STIR ? (ne + 10* prj.NbOS)/104857.6 : ne/131072), 2); //... memory allocation for wm float arrays ................................... - for( int i = 0 ; i < wmh.prj.NbOS ; i++ ){ + for( int i = 0 ; i < this->wmh.prj.NbOS ; i++ ){ - if ( ( wm.val[ i ] = new (std::nothrow) float [ NITEMS[kOS][ i ] ]) == NULL) + if ( ( this->wm.val[ i ] = new (std::nothrow) float [ NITEMS[kOS][ i ] ]) == NULL) { //error_wm_SPECT( 200, "wm.val[][]" ); error("Error allocating space to store values for SPECTUB matrix"); } - if ( ( wm.col[ i ] = new (std::nothrow) int [ NITEMS[kOS][ i ] ]) == NULL) + if ( ( this->wm.col[ i ] = new (std::nothrow) int [ NITEMS[kOS][ i ] ]) == NULL) { //error_wm_SPECT( 200, "wm.col[]" ); error("Error allocating space to store column indices for SPECTUB matrix"); @@ -799,47 +793,47 @@ compute_one_subset(const int kOS, //... to initialize wm to zero ...................... - for ( int i = 0 ; i < wm.NbOS ; i++ ){ + for ( int i = 0 ; i < this->wm.NbOS ; i++ ){ - wm.ne[ i ] = 0; + this->wm.ne[ i ] = 0; for( int j = 0 ; j < NITEMS[kOS][ i ] ; j++ ){ - wm.val[ i ][ j ] = (float)0.; - wm.col[ i ][ j ] = 0; + this->wm.val[ i ][ j ] = (float)0.; + this->wm.col[ i ][ j ] = 0; } } - wm.ne[ wm.NbOS ] = 0; + this->wm.ne[ this->wm.NbOS ] = 0; //... wm calculation for this subset ........................... wm_calculation ( kOS, ang, vox, bin, vol, prj, attmap, msk_3d, msk_2d, maxszb, &gaussdens, NITEMS[kOS], - wm, wmh, Rrad ); + this->wm, this->wmh, Rrad ); info(boost::format("Weight matrix calculation done. time %1% (s)") % timer.value(), 2); //... fill lor ......................... - for( int j = 0 ; j < wm.NbOS ; j++ ){ + for( int j = 0 ; j < this->wm.NbOS ; j++ ){ ProjMatrixElemsForOneBin lor; Bin bin; bin.segment_num()=0; - bin.view_num()=wm.na [ j ]; - bin.axial_pos_num()=wm.ns [ j ]; - bin.tangential_pos_num()=wm.nb [ j ]; + bin.view_num()=this->wm.na [ j ]; + bin.axial_pos_num()=this->wm.ns [ j ]; + bin.tangential_pos_num()=this->wm.nb [ j ]; bin.set_bin_value(0); lor.set_bin(bin); - lor.reserve(wm.ne[ j ]); - for ( int i = 0 ; i < wm.ne[ j ] ; i++ ){ + lor.reserve(this->wm.ne[ j ]); + for ( int i = 0 ; i < this->wm.ne[ j ] ; i++ ){ const ProjMatrixElemsForOneBin::value_type - elem(Coordinate3D(wm.nz[ wm.col[ j ][ i ] ],wm.ny[ wm.col[ j ][ i ] ],wm.nx[ wm.col[ j ][ i ] ]), wm.val[ j ][ i ]); + elem(Coordinate3D(this->wm.nz[ this->wm.col[ j ][ i ] ],this->wm.ny[ this->wm.col[ j ][ i ] ],this->wm.nx[ this->wm.col[ j ][ i ] ]), this->wm.val[ j ][ i ]); lor.push_back( elem); } - delete [] wm.val[ j ]; - delete [] wm.col[ j ]; + delete [] this->wm.val[ j ]; + delete [] this->wm.col[ j ]; this->cache_proj_matrix_elems_for_one_bin(lor); } @@ -876,7 +870,7 @@ calculate_proj_matrix_elems_for_one_bin(ProjMatrixElemsForOneBin& lor } info(boost::format("Computing matrix elements for view %1%") % view_num, 2);// potentially pass a wm, wmh[threadh] not sure if works then in setup we need an array of wmh - compute_one_subset(kOS,wm,wmh,Rrad); + compute_one_subset(kOS,Rrad); subset_already_processed[kOS]=true; } lor.erase(); diff --git a/src/recon_buildblock/SPECTUB_Weight3d.cxx b/src/recon_buildblock/SPECTUB_Weight3d.cxx index aeb6195d5d..06218b9214 100644 --- a/src/recon_buildblock/SPECTUB_Weight3d.cxx +++ b/src/recon_buildblock/SPECTUB_Weight3d.cxx @@ -46,16 +46,16 @@ using namespace std; //========================================================================== void wm_calculation(const int kOS, - const angle_type *const ang, + const angle_type * const ang, voxel_type vox, bin_type bin, - const volume_type& vol, + const volume_type &vol, const proj_type& prj, const float *attmap, const bool *msk_3d, const bool *msk_2d, const int maxszb, - const discrf_type *const gaussdens, + const discrf_type * const gaussdens, const int *const NITEMS, wm_da_type &wm, wmh_type &wmh, float *Rrad) { @@ -174,7 +174,7 @@ void wm_calculation(const int kOS, //... to project voxels onto the detection plane and to calculate other distances ..... - voxel_projection( &vox , &eff , prj.lngcmd2 ); + voxel_projection( &vox , &eff , prj.lngcmd2, wmh ); //... setting PSF to zero ......................................... @@ -186,13 +186,13 @@ void wm_calculation(const int kOS, //... correction for PSF .............................. - if ( !wmh.do_psf ) fill_psf_no ( &psf, &psf1d_h, vox, &ang[ ka ], bin.szdx ); + if ( !wmh.do_psf ) fill_psf_no ( &psf, &psf1d_h, vox, &ang[ ka ], bin.szdx, wmh); else{ - if ( wmh.do_psf_3d ) fill_psf_3d ( &psf, &psf1d_h, &psf1d_v, vox, gaussdens, bin.szdx, bin.thdx, bin.thcmd2 ); + if ( wmh.do_psf_3d ) fill_psf_3d ( &psf, &psf1d_h, &psf1d_v, vox, gaussdens, bin.szdx, bin.thdx, bin.thcmd2, wmh ); - else fill_psf_2d ( &psf, &psf1d_h, vox, gaussdens, bin.szdx ); + else fill_psf_2d ( &psf, &psf1d_h, vox, gaussdens, bin.szdx, wmh); } //... correction for attenuation ................................................. @@ -382,7 +382,7 @@ void wm_size_estimation (int kOS, //... to project voxels onto the detection plane and to calculate other distances ..... - voxel_projection( &vox , &eff , prj.lngcmd2 ); + voxel_projection( &vox , &eff , prj.lngcmd2, wmh); //... setting PSF to zero ......................................... @@ -393,13 +393,13 @@ void wm_size_estimation (int kOS, //... correction for PSF .............................. - if ( !wmh.do_psf ) fill_psf_no ( &psf, &psf1d_h, vox, &ang[ ka ], bin.szdx ); + if ( !wmh.do_psf ) fill_psf_no ( &psf, &psf1d_h, vox, &ang[ ka ], bin.szdx, wmh); else{ - if ( wmh.do_psf_3d ) fill_psf_3d ( &psf, &psf1d_h, &psf1d_v, vox, gaussdens, bin.szdx, bin.thdx, bin.thcmd2 ); + if ( wmh.do_psf_3d ) fill_psf_3d ( &psf, &psf1d_h, &psf1d_v, vox, gaussdens, bin.szdx, bin.thdx, bin.thcmd2, wmh ); - else fill_psf_2d ( &psf, &psf1d_h, vox, gaussdens, bin.szdx ); + else fill_psf_2d ( &psf, &psf1d_h, vox, gaussdens, bin.szdx, wmh); } @@ -528,7 +528,7 @@ void calc_vxprj( angle_type *ang ) //=== voxel_projection ===================================================== //========================================================================== -void voxel_projection ( voxel_type *vox, float * eff, float lngcmd2) +void voxel_projection ( voxel_type *vox, float * eff, float lngcmd2, wmh_type &wmh) { if ( wmh.COL.do_fb ){ // fan_beam @@ -559,7 +559,12 @@ void voxel_projection ( voxel_type *vox, float * eff, float lngcmd2) //=== fill_psf_no ========================================================== //========================================================================== -void fill_psf_no( psf2da_type *psf, psf1d_type * psf1d_h, const voxel_type& vox, angle_type const *const ang , float szdx ) +void fill_psf_no( psf2da_type *psf, + psf1d_type * psf1d_h, + const voxel_type& vox, + angle_type const *const ang , + float szdx, + wmh_type &wmh) { psf1d_h->sgmcm = vox.szcm; @@ -573,7 +578,7 @@ void fill_psf_no( psf2da_type *psf, psf1d_type * psf1d_h, const voxel_type& vox, psf1d_h->lngcmd2 = psf1d_h->lngcm / (float)2.; psf1d_h->efres = ang->vxprj.res * psf1d_h->sgmcm; // to resize discretization resolution once applied sgmcm - calc_psf_bin( vox.xd0, wmh.prj.szcm, &ang->vxprj, psf1d_h ); + calc_psf_bin( vox.xd0, wmh.prj.szcm, &ang->vxprj, psf1d_h, wmh); for ( int ie = 0 ; ie < psf1d_h->Nib ; ie++ ){ @@ -588,7 +593,7 @@ void fill_psf_no( psf2da_type *psf, psf1d_type * psf1d_h, const voxel_type& vox, //=== fill_psf_2d ========================================================== //========================================================================== -void fill_psf_2d( psf2da_type *psf, psf1d_type * psf1d_h, const voxel_type& vox, discrf_type const* const gaussdens, float szdx ) +void fill_psf_2d( psf2da_type *psf, psf1d_type * psf1d_h, const voxel_type& vox, discrf_type const* const gaussdens, float szdx,wmh_type &wmh ) { psf1d_h->sgmcm = calc_sigma_h( vox, wmh.COL ); @@ -599,7 +604,7 @@ void fill_psf_2d( psf2da_type *psf, psf1d_type * psf1d_h, const voxel_type& vox, psf1d_h->efres = gaussdens->res * psf1d_h->sgmcm ; - calc_psf_bin( vox.xd0, wmh.prj.szcm, gaussdens, psf1d_h ); + calc_psf_bin( vox.xd0, wmh.prj.szcm, gaussdens, psf1d_h, wmh); for ( int ie = 0 ; ie < psf1d_h->Nib ; ie++ ){ @@ -618,7 +623,10 @@ void fill_psf_2d( psf2da_type *psf, psf1d_type * psf1d_h, const voxel_type& vox, void fill_psf_3d (psf2da_type *psf, psf1d_type *psf1d_h, psf1d_type *psf1d_v, - const voxel_type& vox, discrf_type const * const gaussdens, float szdx, float thdx, float thcmd2 ) + const voxel_type& vox, + discrf_type const * const gaussdens, + float szdx, float thdx, float thcmd2, + SPECTUB::wmh_type &wmh) { //... horizontal component ........................... @@ -638,7 +646,7 @@ void fill_psf_3d (psf2da_type *psf, //... calculation of the horizontal component of psf ................... - calc_psf_bin( vox.xd0, wmh.prj.szcm, gaussdens, psf1d_h ); + calc_psf_bin( vox.xd0, wmh.prj.szcm, gaussdens, psf1d_h, wmh); //... vertical component .............................. @@ -657,7 +665,7 @@ void fill_psf_3d (psf2da_type *psf, //... calculation of the vertical component of psf .................... - calc_psf_bin( thcmd2, wmh.prj.thcm, gaussdens, psf1d_v ); + calc_psf_bin( thcmd2, wmh.prj.thcm, gaussdens, psf1d_v, wmh); //... mixing and setting PSF area to 1 (to correct for tail truncation of Gaussian function) ..... @@ -701,7 +709,8 @@ void fill_psf_3d (psf2da_type *psf, void calc_psf_bin (float center_psf, float binszcm, discrf_type const * const vxprj, - psf1d_type *psf) + psf1d_type *psf, + SPECTUB::wmh_type &wmh) { float weight, preval; From 573e4bda60e20fadc65289405350b3e106465748 Mon Sep 17 00:00:00 2001 From: nmdicom-recon Date: Tue, 18 Apr 2023 15:15:45 +0100 Subject: [PATCH 3/4] remove dependence from global variable; use const * Rrad where possible. --- .../recon_buildblock/ProjMatrixByBinSPECTUB.h | 2 +- .../stir/recon_buildblock/SPECTUB_Tools.h | 22 ++++++------- .../stir/recon_buildblock/SPECTUB_Weight3d.h | 20 ++++++------ .../ProjMatrixByBinSPECTUB.cxx | 12 +++---- src/recon_buildblock/SPECTUB_Tools.cxx | 32 +++++++++---------- src/recon_buildblock/SPECTUB_Weight3d.cxx | 6 ++-- 6 files changed, 47 insertions(+), 47 deletions(-) diff --git a/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h b/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h index 4127343d01..3ca43a9683 100644 --- a/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h +++ b/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h @@ -212,7 +212,7 @@ class ProjMatrixByBinSPECTUB : void compute_one_subset(const int kOS, - float * Rrad) const; + const float *Rrad) const; void delete_UB_SPECT_arrays(); mutable std::vector subset_already_processed; }; diff --git a/src/include/stir/recon_buildblock/SPECTUB_Tools.h b/src/include/stir/recon_buildblock/SPECTUB_Tools.h index 2f7d738ab0..467aae3e36 100644 --- a/src/include/stir/recon_buildblock/SPECTUB_Tools.h +++ b/src/include/stir/recon_buildblock/SPECTUB_Tools.h @@ -319,16 +319,16 @@ typedef struct //... functions from wmtools_SPECT.cpp ......................................... -void write_wm_FC (); // to write double array weight matrix +void write_wm_FC (wm_da_type &wm); // to write double array weight matrix -void write_wm_hdr (); // to write header of a matrix +void write_wm_hdr (wm_da_type &wm, wmh_type &wmh); // to write header of a matrix -void write_wm_STIR (); // to write matrix in STIR format +void write_wm_STIR (SPECTUB::wm_da_type& wm); // to write matrix in STIR format -void index_calc ( int *indexs ); // to calculate projection index order in subsets +void index_calc (int *indexs , wmh_type &wmh); // to calculate projection index order in subsets -void read_Rrad (); // to read variable rotation radius from a text file (1 radius per line) +void read_Rrad (float *Rrad, wmh_type &wmh); // to read variable rotation radius from a text file (1 radius per line) // //void col_params ( collim_type *COL ); // to fill collimator structure @@ -336,21 +336,21 @@ void read_Rrad (); // to read variable rotation radius from a text fil //void read_col_params ( collim_type *COL); // to read collimator parameters from a file -void fill_ang ( angle_type *ang ); // to fill angle structure +void fill_ang (angle_type *ang , wmh_type &wmh, const float *Rrad); // to fill angle structure -void generate_msk ( bool *msk_3d, bool *msk_2d, float *att, volume_type *vol); // to create a boolean mask for wm (no weights outside the msk) +void generate_msk (bool *msk_3d, bool *msk_2d, float *att, volume_type *vol, wmh_type &wmh); // to create a boolean mask for wm (no weights outside the msk) -void read_msk_file ( bool * msk ); // to read mask from a file +void read_msk_file (bool * msk , wmh_type &wmh); // to read mask from a file -void read_att_map ( float *attmap ); // to read attenuation map from a file +void read_att_map (float *attmap , wmh_type &wmh); // to read attenuation map from a file -int max_psf_szb ( angle_type *ang ); +int max_psf_szb (angle_type *ang , wmh_type &wmh); float calc_sigma_h ( voxel_type vox, collim_type COL); -float calc_sigma_v ( voxel_type vox, collim_type COL); +float calc_sigma_v (voxel_type vox, collim_type COL, wmh_type &wmh); char *itoa ( int n, char *s); // to conver integer to ascii diff --git a/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h b/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h index 57888f662f..5ac3f71cb8 100644 --- a/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h +++ b/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h @@ -32,22 +32,22 @@ namespace SPECTUB { // extern float * Rrad; //! radii per view -void wm_calculation( const int kOS, +void wm_calculation(const int kOS, const SPECTUB::angle_type *const ang, SPECTUB::voxel_type vox, - bin_type bin, + bin_type bin, const SPECTUB::volume_type& vol, - const proj_type& prj, - const float *attmap, - const bool *msk_3d, - const bool *msk_2d, - const int maxszb, + const proj_type& prj, + const float *attmap, + const bool *msk_3d, + const bool *msk_2d, + const int maxszb, const SPECTUB::discrf_type *const gaussdens, const int *const NITEMS, SPECTUB::wm_da_type& wm, SPECTUB::wmh_type& wmh, - float * Rrad - ); + const float *Rrad + ); void wm_size_estimation (int kOS, const SPECTUB::angle_type * const ang, @@ -61,7 +61,7 @@ void wm_size_estimation (int kOS, const SPECTUB::discrf_type * const gaussdens, int *NITEMS, SPECTUB::wmh_type &wmh, - float *Rrad); + const float *Rrad); //... geometric component ............................................ diff --git a/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx b/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx index 79f9066a56..7c5d263c2d 100644 --- a/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx +++ b/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx @@ -502,12 +502,12 @@ set_up( //... to sort angles into subsets ...................................... prj.order = new int [ prj.Nang ]; - index_calc( prj.order ); + index_calc( prj.order, wmh); //... to fill ang structure ............................................ ang = new angle_type [ prj.Nang ]; - fill_ang( ang ); + fill_ang( ang,wmh, Rrad ); //... to fill high resolution discrete distribution functions .............. @@ -574,12 +574,12 @@ set_up( // we do this to avoid using its own read_msk_file wmh.do_msk_file = false; wmh.do_msk_att = true; - generate_msk( msk_3d, msk_2d, mask_from_file, &vol); + generate_msk( msk_3d, msk_2d, mask_from_file, &vol, wmh); delete[] mask_from_file; } else { - generate_msk( msk_3d, msk_2d, attmap, &vol); + generate_msk( msk_3d, msk_2d, attmap, &vol, wmh); } } else msk_2d = msk_3d = NULL; @@ -591,7 +591,7 @@ set_up( //... setting PSF maximum size (in bins) and memory allocation for PSF values ....... - this->maxszb = max_psf_szb( ang ); // maximum PSF size (horizontal component of PSF) + this->maxszb = max_psf_szb( ang, wmh ); // maximum PSF size (horizontal component of PSF) NITEMS = new int * [prj.NOS]; for (int kOS=0; kOS() / (float)180. ; @@ -614,7 +614,7 @@ void fill_ang ( angle_type *ang ) //=== generate msk ============================================================ //============================================================================= -void generate_msk ( bool *msk_3d, bool *msk_2d, float *attmap, volume_type * vol ) +void generate_msk ( bool *msk_3d, bool *msk_2d, float *attmap, volume_type * vol,SPECTUB::wmh_type& wmh ) { //... initialzation of msk to true ......................... @@ -672,7 +672,7 @@ void generate_msk ( bool *msk_3d, bool *msk_2d, float *attmap, volume_type * vol else { //... to read a mask from a (int) file .................... - if ( wmh.do_msk_file ) read_msk_file( msk_3d ); + if ( wmh.do_msk_file ) read_msk_file( msk_3d, wmh); } } @@ -717,7 +717,7 @@ void generate_msk ( bool *msk_3d, bool *msk_2d, float *attmap, volume_type * vol //=== read_mask file ========================================================== //============================================================================= -void read_msk_file( bool *msk ) +void read_msk_file( bool *msk, SPECTUB::wmh_type& wmh ) { FILE *fid; int *aux; @@ -739,7 +739,7 @@ void read_msk_file( bool *msk ) //=== read_att_map ============================================================ //============================================================================= -void read_att_map( float *attmap ) +void read_att_map( float *attmap,SPECTUB::wmh_type& wmh ) { FILE *fid; if ( ( fid = fopen( wmh.att_fn.c_str() , "rb") ) == NULL ) error_wmtools_SPECT ( 124, wmh.att_fn ); @@ -763,7 +763,7 @@ void read_att_map( float *attmap ) //=== max_psf_szb ========================================================== //========================================================================== -int max_psf_szb( angle_type *ang ) +int max_psf_szb( angle_type *ang, SPECTUB::wmh_type& wmh) { int maxszb; float Rrad_max = ang[0].Rrad; @@ -803,7 +803,7 @@ int max_psf_szb( angle_type *ang ) maxszb = (int)floor( wmh.maxsigm * (float)2. * sig_h_max_cm / wmh.prj.szcm ) + 3; if ( wmh.do_psf_3d ){ - float sig_v_max_cm = calc_sigma_v( vox, wmh.COL ); + float sig_v_max_cm = calc_sigma_v( vox, wmh.COL, wmh ); int maxszb_v = (int)floor( wmh.maxsigm * (float)2. * sig_v_max_cm / wmh.prj.thcm ) + 3; maxszb = maxim( maxszb , maxszb_v ); } diff --git a/src/recon_buildblock/SPECTUB_Weight3d.cxx b/src/recon_buildblock/SPECTUB_Weight3d.cxx index 06218b9214..c254004a5b 100644 --- a/src/recon_buildblock/SPECTUB_Weight3d.cxx +++ b/src/recon_buildblock/SPECTUB_Weight3d.cxx @@ -56,7 +56,7 @@ void wm_calculation(const int kOS, const bool *msk_2d, const int maxszb, const discrf_type * const gaussdens, - const int *const NITEMS, wm_da_type &wm, wmh_type &wmh, float *Rrad) + const int *const NITEMS, wm_da_type &wm, wmh_type &wmh, const float *Rrad) { float weight; @@ -314,7 +314,7 @@ void wm_size_estimation (int kOS, const discrf_type * const gaussdens, int *NITEMS, wmh_type& wmh, - float * Rrad) + const float * Rrad) { int jp; float eff; @@ -650,7 +650,7 @@ void fill_psf_3d (psf2da_type *psf, //... vertical component .............................. - psf1d_v->sgmcm = calc_sigma_v( vox, wmh.COL); + psf1d_v->sgmcm = calc_sigma_v( vox, wmh.COL, wmh); psf1d_v->lngcmd2 = psf1d_v->sgmcm * wmh.maxsigm; psf1d_v->lngcm = psf1d_v->lngcmd2 * (float)2.; psf1d_v->di = min( (int) floor( thdx / psf1d_v->sgmcm ), gaussdens->lng - 1 ) ; From 581c302f8232e161f09a6b83391939816795f530 Mon Sep 17 00:00:00 2001 From: nmdicom-recon Date: Wed, 19 Apr 2023 10:21:49 +0100 Subject: [PATCH 4/4] remove comments and add info to release5.3 --- documentation/release_5.2.htm | 3 +++ src/include/stir/recon_buildblock/SPECTUB_Weight3d.h | 5 ----- src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx | 11 ----------- src/recon_buildblock/SPECTUB_Tools.cxx | 8 -------- 4 files changed, 3 insertions(+), 24 deletions(-) diff --git a/documentation/release_5.2.htm b/documentation/release_5.2.htm index 82ac1762cc..ba6d283e49 100644 --- a/documentation/release_5.2.htm +++ b/documentation/release_5.2.htm @@ -31,6 +31,9 @@

Bug fixes

New functionality

    +
  • Global variables in SPECTUB have been substituted by class members, such that multiple SPECTUB projectors can be used. +
    PR #1169. +
  • Scatter estimation is now smoothed in axial direction for BlocksOnCylindrical scanners.
    PR #1172.
    • diff --git a/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h b/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h index 5ac3f71cb8..8f15d93467 100644 --- a/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h +++ b/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h @@ -27,11 +27,6 @@ namespace SPECTUB { -// extern wm_da_type wm; //! weight (or probability) matrix -// extern wmh_type wmh; //! information to construct wm -// extern float * Rrad; //! radii per view - - void wm_calculation(const int kOS, const SPECTUB::angle_type *const ang, SPECTUB::voxel_type vox, diff --git a/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx b/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx index 7c5d263c2d..8516d6de89 100644 --- a/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx +++ b/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx @@ -21,14 +21,11 @@ \author Daniel Deidda */ -//#include "stir/ProjDataInterfile.h" #include "stir/recon_buildblock/ProjMatrixByBinSPECTUB.h" #include "stir/recon_buildblock/TrivialDataSymmetriesForBins.h" #include "stir/ProjDataInfoCylindricalArcCorr.h" -//#include "stir/KeyParser.h" #include "stir/IO/read_from_file.h" #include "stir/ProjDataInfo.h" -//#include "stir/utilities.h" #include "stir/VoxelsOnCartesianGrid.h" #include "stir/Succeeded.h" #include "stir/is_null_ptr.h" @@ -55,19 +52,11 @@ #include #include #include -//#include //... user defined libraries ............................................................. #include "stir/recon_buildblock/SPECTUB_Weight3d.h" -/* UB-SPECT global variables */ -//namespace SPECTUB { -// wm_da_type wm; -// wmh_type wmh; -// float * Rrad; -//} - using namespace SPECTUB; START_NAMESPACE_STIR diff --git a/src/recon_buildblock/SPECTUB_Tools.cxx b/src/recon_buildblock/SPECTUB_Tools.cxx index 33d699955a..e71a9f0e70 100644 --- a/src/recon_buildblock/SPECTUB_Tools.cxx +++ b/src/recon_buildblock/SPECTUB_Tools.cxx @@ -22,7 +22,6 @@ #include using namespace std; -//using std::string; //... user defined libraries ....................................... @@ -45,13 +44,6 @@ namespace SPECTUB { #define DELIMITER1 '#' //delimiter character in input parameter text file #define DELIMITER2 '%' //delimiter character in input parameter text file -//... global variables .............................................. - -//extern wm_da_type wm; -//extern wmh_type wmh; -//extern float * Rrad; - - //============================================================================= //=== write_wm_FC ============================================================= //=============================================================================