Skip to content

Commit

Permalink
fix SPECT subsets functionality - not perfect yet - tests need to be …
Browse files Browse the repository at this point in the history 
…done properly
  • Loading branch information
Sam Porter committed Nov 20, 2024
1 parent d74fcc1 commit 3bf6307
Show file tree
Hide file tree
Showing 5 changed files with 870 additions and 186 deletions.
7 changes: 3 additions & 4 deletions src/include/stir/recon_buildblock/SPECTUB_Tools.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,10 +87,9 @@ typedef struct
int Nsli; // number of slices
float thcm; // slice thickness in cm

int Nang; // number of projection angles
float ang0; // initial projection angle. degrees from upper detection plane (parallel to table). Negative for CW rotacions (see
// manual)
float incr; // angle increment between two consecutive projection angles. Degrees. Negative for CW, Positive for CCW
int Nang; // number of projection angles
std::vector<float> angles; // projection angles. degrees from upper detection plane (parallel to table). Negative for CW
// rotations (see manual)

int NOS; // number of subsets in which to split the matrix
int NangOS; // Number of angles in each subset = Nang/NOS
Expand Down
54 changes: 44 additions & 10 deletions src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,12 +24,14 @@
#include "stir/recon_buildblock/ProjMatrixByBinSPECTUB.h"
#include "stir/recon_buildblock/TrivialDataSymmetriesForBins.h"
#include "stir/ProjDataInfoCylindricalArcCorr.h"
#include "stir/ProjDataInfoSubsetByView.h"
#include "stir/IO/read_from_file.h"
#include "stir/ProjDataInfo.h"
#include "stir/VoxelsOnCartesianGrid.h"
#include "stir/Succeeded.h"
#include "stir/is_null_ptr.h"
#include "stir/Coordinate3D.h"
#include "stir/stream.h"
#include "stir/info.h"
#include "stir/warning.h"
#include "stir/error.h"
Expand Down Expand Up @@ -256,9 +258,6 @@ ProjMatrixByBinSPECTUB::set_up(const shared_ptr<const ProjDataInfo>& proj_data_i
this->voxel_size = image_info_ptr->get_voxel_size();
this->origin = image_info_ptr->get_origin();

const ProjDataInfoCylindricalArcCorr* proj_Data_Info_Cylindrical
= dynamic_cast<const ProjDataInfoCylindricalArcCorr*>(this->proj_data_info_ptr.get());

CPUTimer timer;
timer.start();

Expand Down Expand Up @@ -301,12 +300,20 @@ ProjMatrixByBinSPECTUB::set_up(const shared_ptr<const ProjDataInfo>& proj_data_i
vox.thcm = vol.thcm;

//... projecction parameters ..........................................
prj.ang0 = this->proj_data_info_ptr->get_scanner_ptr()->get_intrinsic_azimuthal_tilt() * float(180 / _PI);
prj.incr = proj_Data_Info_Cylindrical->get_azimuthal_angle_sampling() * float(180 / _PI);
prj.thcm = proj_Data_Info_Cylindrical->get_axial_sampling(0) / 10;
prj.angles = std::vector<float>(prj.Nang);
{
Bin bin;
for (int i = 0; i < prj.Nang; i++)
{
bin.view_num() = i;
prj.angles[i] = static_cast<float>(this->proj_data_info_ptr->get_phi(bin) * 180 / _PI);
}
// all bins will have same axial sampling
prj.thcm = this->proj_data_info_ptr->get_sampling_in_m(bin) / 10;
}

//.......geometrical and other derived parameters of projection structure...........
prj.Nsli = proj_Data_Info_Cylindrical->get_num_axial_poss(0); // number of slices
prj.Nsli = this->proj_data_info_ptr->get_num_axial_poss(0); // number of slices
prj.lngcm = prj.Nbin * prj.szcm; // length in cm of the detection line
prj.Nbp = prj.Nbin * prj.Nsli; // number of bins for each projection angle (2D-projection)
prj.Nbt = prj.Nbp * prj.Nang; // total number of bins considering all the projection angles
Expand All @@ -330,7 +337,35 @@ ProjMatrixByBinSPECTUB::set_up(const shared_ptr<const ProjDataInfo>& proj_data_i
"SPECTUB Matrix (probably) only works with equal number of slices for projection data (%1%) and image (%2%)")
% wmh.prj.Nsli % vol.Nsli);
//....rotation radius .................................................
const VectorWithOffset<float> radius_all_views = proj_Data_Info_Cylindrical->get_ring_radii_for_all_views();
VectorWithOffset<float> radius_all_views(prj.Nang);
{
if (auto proj_Data_Info_Cylindrical = dynamic_cast<const ProjDataInfoCylindricalArcCorr*>(this->proj_data_info_ptr.get()))
{
radius_all_views = proj_Data_Info_Cylindrical->get_ring_radii_for_all_views();
}
else if (auto proj_data_info_subset_ptr = dynamic_cast<const ProjDataInfoSubsetByView *>(this->proj_data_info_ptr.get()))
{
if (auto proj_Data_Info_Cylindrical = dynamic_cast<const ProjDataInfoCylindricalArcCorr*>(
proj_data_info_subset_ptr->get_original_proj_data_info_sptr().get()))
{
for (int i = 0; i < prj.Nang; ++i)
{
Bin bin;
bin.view_num() = i;
const int org_view = proj_data_info_subset_ptr->get_original_bin(bin).view_num();
radius_all_views[i] = proj_Data_Info_Cylindrical->get_ring_radius(org_view);
}
}
else
{
error("ProjMatrixByBinSPECTUB: can only handle ProjDataInfoCylindricalArcCorr");
}
}
else
{
error("ProjMatrixByBinSPECTUB: can only handle ProjDataInfoCylindricalArcCorr");
}
}

{
const auto max_radius = *std::max_element(radius_all_views.begin(), radius_all_views.end());
Expand Down Expand Up @@ -494,8 +529,7 @@ ProjMatrixByBinSPECTUB::set_up(const shared_ptr<const ProjDataInfo>& proj_data_i
info_stream << "Number of slices: " << wmh.vol.Nsli << "\tslice_thickness: " << wmh.vol.thcm << std::endl;
info_stream << "Number of bins: " << wmh.prj.Nbin << "\tbin size: " << wmh.prj.szcm << "\taxial size: " << wmh.prj.thcm
<< std::endl;
info_stream << "Number of angles: " << wmh.prj.Nang << "\tAngle increment: " << wmh.prj.incr
<< "\tFirst angle: " << wmh.prj.ang0 << std::endl;
info_stream << "Number of angles: " << wmh.prj.Nang << "\tangles (deg): " << wmh.prj.angles << std::endl;
info_stream << "Number of subsets: " << wmh.prj.NOS << std::endl;
if (wmh.do_att)
{
Expand Down
10 changes: 7 additions & 3 deletions src/recon_buildblock/SPECTUB_Tools.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include <math.h>
#include <boost/math/special_functions/fpclassify.hpp>
#include "stir/error.h"
#include "stir/stream.h"
#include <boost/format.hpp>
#include <boost/math/constants/constants.hpp>

Expand Down Expand Up @@ -129,9 +130,12 @@ write_wm_hdr(SPECTUB::wm_da_type& wm, SPECTUB::wmh_type& wmh)
stream1 << "number of bins per line: " << wmh.prj.Nbin << endl;
stream1 << "bin size (cm): " << wmh.prj.szcm << endl;
stream1 << "number of angles: " << wmh.prj.Nang << endl;
stream1 << "first angle (deg): " << wmh.prj.ang0 << endl;
stream1 << "angle increment between consecutive projections (deg): " << wmh.prj.incr << endl;

{
using namespace stir;
stream1 << "angles (deg): " << wmh.prj.angles << endl;
}

stream1 << "first slice to reconstruct : " << wmh.vol.first_sl << endl;
stream1 << "last slice to reconstruct : " << wmh.vol.last_sl << endl;
stream1 << "number of subsets in which to split the matrix: " << wmh.prj.NOS << endl;
Expand Down Expand Up @@ -606,7 +610,7 @@ fill_ang(angle_type* ang, SPECTUB::wmh_type& wmh, const float* Rrad)

//... ratios calculation .......................................................

float deg = wmh.prj.ang0 + (float)i * wmh.prj.incr; // angle in degrees
const float deg = wmh.prj.angles[i];
ang[i].cos = cos(deg * dg2rd); // cosinus of the angle
ang[i].sin = sin(deg * dg2rd); // sinus of the angle

Expand Down
Loading

0 comments on commit 3bf6307

Please sign in to comment.