Skip to content

Commit

Permalink
Added an additional test that compares sinogram bins for LORs that do…
Browse files Browse the repository at this point in the history
… not change during scanner downsampling.
  • Loading branch information
Markus Jehl committed Aug 30, 2024
1 parent ee436f9 commit f3e697d
Show file tree
Hide file tree
Showing 2 changed files with 155 additions and 3 deletions.
2 changes: 1 addition & 1 deletion src/test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ set(buildblock_simple_tests
include(stir_test_exe_targets)

foreach(source ${buildblock_simple_tests})
create_stir_test(${source} "buildblock;IO;buildblock;numerics_buildblock;display;IO;recon_buildblock;Shape_buildblock" "")
create_stir_test(${source} "buildblock;IO;buildblock;numerics_buildblock;display;IO;recon_buildblock;Shape_buildblock;scatter_buildblock" "")
endforeach()
#

Expand Down
156 changes: 154 additions & 2 deletions src/test/test_interpolate_projdata.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@
#include "stir/Shape/Box3D.h"
#include "stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h"
#include "stir/recon_buildblock/ForwardProjectorByBinUsingProjMatrixByBin.h"
#include "stir/scatter/SingleScatterSimulation.h"

#include "stir/RunTests.h"

Expand All @@ -57,6 +58,7 @@ class InterpolationTests : public RunTests
void scatter_interpolation_test_blocks_asymmetric();
void scatter_interpolation_test_cyl_asymmetric();
void scatter_interpolation_test_blocks_downsampled();
void transaxial_upsampling_interpolation_test_blocks();

void check_symmetry(const SegmentBySinogram<float>& segment);
void compare_segment(const SegmentBySinogram<float>& segment1, const SegmentBySinogram<float>& segment2, float maxDiff);
Expand Down Expand Up @@ -754,7 +756,7 @@ InterpolationTests::scatter_interpolation_test_blocks_downsampled()
auto proj_data = ProjDataInMemory(std::make_shared<ExamInfo>(exam_info), proj_data_info);
auto downsampled_proj_data = ProjDataInMemory(std::make_shared<ExamInfo>(exam_info), downsampled_proj_data_info);

// define a cylinder precisely in the middle of the FOV, such that symmetry can be used for validation
// define a cylinder and a box that are off-centre, such that the shapes in the sinogram can be compared
auto emission_map = VoxelsOnCartesianGrid<float>(*proj_data_info, 1);
auto cyl_map = VoxelsOnCartesianGrid<float>(*proj_data_info, 1);
auto cylinder = EllipsoidalCylinder(40, 40, 20, CartesianCoordinate3D<float>(80, 100, 0));
Expand All @@ -763,7 +765,7 @@ InterpolationTests::scatter_interpolation_test_blocks_downsampled()
box.construct_volume(emission_map, CartesianCoordinate3D<int>(1, 1, 1));
emission_map += cyl_map;

// project the cylinder onto the full-scale scanner proj data
// project the emission map onto the full-scale scanner proj data
auto pm = ProjMatrixByBinUsingRayTracing();
pm.set_use_actual_detector_boundaries(true);
pm.enable_cache(false);
Expand Down Expand Up @@ -801,6 +803,155 @@ InterpolationTests::scatter_interpolation_test_blocks_downsampled()
compare_segment_shape(full_size_model_sino.get_segment_by_sinogram(0), interpolated_proj_data.get_segment_by_sinogram(0), 3);
}

void
InterpolationTests::transaxial_upsampling_interpolation_test_blocks()
{
info("Performing transaxial downampled interpolation test for BlocksOnCylindrical scanner");
auto time_frame_def = TimeFrameDefinitions();
time_frame_def.set_num_time_frames(1);
time_frame_def.set_time_frame(1, 0, 1e9);
auto exam_info = ExamInfo();
exam_info.set_high_energy_thres(650);
exam_info.set_low_energy_thres(425);
exam_info.set_time_frame_definitions(time_frame_def);

// define the original scanner and a downsampled one, as it would be used for scatter simulation
auto scanner = Scanner(Scanner::User_defined_scanner,
"Some_BlocksOnCylindrical_Scanner",
96,
3,
60,
60,
127,
6.5,
3.313,
1.65,
-3.1091819,
1,
3,
3,
4,
1,
1,
1,
0.14,
511,
1,
0,
500,
"BlocksOnCylindrical",
3.313,
7.0,
20.0,
29.0);
auto proj_data_info = shared_ptr<ProjDataInfo>(
std::move(ProjDataInfo::construct_proj_data_info(std::make_shared<Scanner>(scanner), 1, 0, 48, 60, false)));

// use the code in scatter simulation to downsample the scanner
auto scatter_simulation = SingleScatterSimulation();
scatter_simulation.set_template_proj_data_info(*proj_data_info);
scatter_simulation.set_exam_info(exam_info);
scatter_simulation.downsample_scanner(-1, 96 / 4); // number of detectors per ring reduced by factor of four
auto downsampled_proj_data_info = scatter_simulation.get_template_proj_data_info_sptr();

auto proj_data = ProjDataInMemory(std::make_shared<ExamInfo>(exam_info), proj_data_info);
auto downsampled_proj_data = ProjDataInMemory(std::make_shared<ExamInfo>(exam_info), downsampled_proj_data_info);

// define a cylinder precisely in the middle of the FOV
auto emission_map = VoxelsOnCartesianGrid<float>(*downsampled_proj_data_info, 1);
make_symmetric_object(emission_map);

// project the cylinder onto the full-scale scanner proj data
auto pm = ProjMatrixByBinUsingRayTracing();
pm.set_use_actual_detector_boundaries(true);
pm.enable_cache(false);
auto forw_proj = ForwardProjectorByBinUsingProjMatrixByBin(std::make_shared<ProjMatrixByBinUsingRayTracing>(pm));
forw_proj.set_up(proj_data_info, std::make_shared<VoxelsOnCartesianGrid<float>>(emission_map));
auto full_size_model_sino = ProjDataInMemory(proj_data);
full_size_model_sino.fill(0);
forw_proj.forward_project(full_size_model_sino, emission_map);

// also project onto the downsampled scanner
emission_map = VoxelsOnCartesianGrid<float>(*downsampled_proj_data_info, 1);
make_symmetric_object(emission_map);
forw_proj.set_up(downsampled_proj_data_info, std::make_shared<VoxelsOnCartesianGrid<float>>(emission_map));
auto downsampled_model_sino = ProjDataInMemory(downsampled_proj_data);
downsampled_model_sino.fill(0);
forw_proj.forward_project(downsampled_model_sino, emission_map);

// write the proj data to file
downsampled_model_sino.write_to_file("transaxially_downsampled_sino_for_LOR.hs");

// interpolate the downsampled proj data to the original scanner size and fill in oblique sinograms
auto interpolated_direct_proj_data = ProjDataInMemory(proj_data);
interpolate_projdata(interpolated_direct_proj_data, downsampled_model_sino, BSpline::linear, false);
auto interpolated_proj_data = ProjDataInMemory(proj_data);
inverse_SSRB(interpolated_proj_data, interpolated_direct_proj_data);

// write the proj data to file
interpolated_proj_data.write_to_file("transaxially_interpolated_sino_for_LOR.hs");

// Identify the bins which should be identical between the downsampled and the interpolated sinogram:
// Each module has 96 / 8 = 12 crystal in the full size scanner, organised in 3 blocks of 4 crystals, while
// the downsampled scanner has 3 crystals per module. The idea is that the centre of the outer two
// is in exactly the same position than the centre of the first and last crystal in the full size scanner.
SegmentBySinogram<float> sinogram_downsampled = downsampled_proj_data.get_empty_segment_by_sinogram(0, false, 0);
SegmentBySinogram<float> sinogram_full_size = proj_data.get_empty_segment_by_sinogram(0, false, 0);
const auto pdi_downsampled
= dynamic_cast<const ProjDataInfoGenericNoArcCorr*>(&(*downsampled_proj_data.get_proj_data_info_sptr()));
const auto pdi_full_size = dynamic_cast<const ProjDataInfoGenericNoArcCorr*>(&(*sinogram_full_size.get_proj_data_info_sptr()));

int tested_LORs = 0;
for (int det1_downsampled = 0; det1_downsampled < 3 * 8; det1_downsampled++)
{
if (det1_downsampled % 3 == 1)
continue; // skip the central crystal of each module
for (int det2_downsampled = 0; det2_downsampled < 3 * 8; det2_downsampled++)
{
if (det2_downsampled % 3 == 1 || det1_downsampled == det2_downsampled)
continue; // skip the central crystal of each module
if (det1_downsampled / 3 == det2_downsampled / 3)
continue; // skip the LORs that lie on the same module

int view_ds, tang_pos_ds;
pdi_downsampled->get_view_tangential_pos_num_for_det_num_pair(view_ds, tang_pos_ds, det1_downsampled, det2_downsampled);
BasicCoordinate<3, int> index_downsampled;
index_downsampled[1] = 1; // looking at central slice
index_downsampled[2] = view_ds;
index_downsampled[3] = tang_pos_ds;

if (tang_pos_ds < pdi_downsampled->get_min_tangential_pos_num()
|| tang_pos_ds > pdi_downsampled->get_max_tangential_pos_num())
continue;

int view_fs, tang_pos_fs;
pdi_full_size->get_view_tangential_pos_num_for_det_num_pair(
view_fs,
tang_pos_fs,
(det1_downsampled / 3) * 12 + ((det1_downsampled % 3) / 2) * 11,
(det2_downsampled / 3) * 12 + ((det2_downsampled % 3) / 2) * 11);

BasicCoordinate<3, int> index_full_size;
index_full_size[1] = 1; // looking at central slice
index_full_size[2] = view_fs;
index_full_size[3] = tang_pos_fs;

if (tang_pos_fs < pdi_full_size->get_min_tangential_pos_num()
|| tang_pos_fs > pdi_full_size->get_max_tangential_pos_num())
continue;

// confirm that the difference is smaller than an empirically found value
check_if_less(std::abs(sinogram_downsampled[index_downsampled] - sinogram_full_size[index_full_size]),
0.01,
"difference between sinogram bin is larger than expected");

tested_LORs++;
}
}

info(boost::format("A total of %1% LORs were compared between the downsampled and the interpolated sinogram.") % tested_LORs);
}

void
InterpolationTests::run_tests()
{
Expand All @@ -809,6 +960,7 @@ InterpolationTests::run_tests()
scatter_interpolation_test_blocks_asymmetric();
scatter_interpolation_test_cyl_asymmetric();
scatter_interpolation_test_blocks_downsampled();
transaxial_upsampling_interpolation_test_blocks();
}

END_NAMESPACE_STIR
Expand Down

0 comments on commit f3e697d

Please sign in to comment.