From 5bedd9321b6e5656d7fba3d840e3cebbc8919c4c Mon Sep 17 00:00:00 2001 From: nmdicom-recon Date: Tue, 3 Oct 2023 13:51:05 +0100 Subject: [PATCH 1/9] changing timing index --- src/buildblock/ProjDataInfo.cxx | 12 ++++++------ src/buildblock/ProjDataInfoCylindrical.cxx | 4 ++-- .../ProjDataInfoCylindricalNoArcCorr.cxx | 16 ++++++++-------- .../stir/ProjDataInfoCylindricalNoArcCorr.inl | 4 ++-- .../DataSymmetriesForBins_PET_CartesianGrid.inl | 12 ++++++++---- .../SymmetryOperations_PET_CartesianGrid.h | 10 ++++++---- .../SymmetryOperations_PET_CartesianGrid.inl | 4 ++-- src/test/test_time_of_flight.cxx | 2 +- 8 files changed, 35 insertions(+), 29 deletions(-) diff --git a/src/buildblock/ProjDataInfo.cxx b/src/buildblock/ProjDataInfo.cxx index f032536f0a..cf391c0447 100644 --- a/src/buildblock/ProjDataInfo.cxx +++ b/src/buildblock/ProjDataInfo.cxx @@ -74,9 +74,9 @@ float ProjDataInfo::get_k(const Bin& bin) const { if (!(num_tof_bins%2)) - return bin.timing_pos_num() * tof_increament_in_mm + tof_increament_in_mm / 2.f; + return (bin.timing_pos_num() - (num_tof_bins)/2.f)* tof_increament_in_mm ; else - return (bin.timing_pos_num() * tof_increament_in_mm); + return (bin.timing_pos_num() - (num_tof_bins-1)/2.f)* tof_increament_in_mm ; } double @@ -228,14 +228,14 @@ ProjDataInfo::set_tof_mash_factor(const int new_num) tof_increament_in_mm = tof_delta_time_to_mm(tof_mash_factor * scanner_ptr->get_size_of_timing_pos()); // TODO cope with even numbers! - min_tof_pos_num = - (scanner_ptr->get_max_num_timing_poss() / tof_mash_factor)/2; - max_tof_pos_num = min_tof_pos_num + (scanner_ptr->get_max_num_timing_poss() / tof_mash_factor) -1; + min_tof_pos_num = 0;//- (scanner_ptr->get_max_num_timing_poss() / tof_mash_factor)/2; + max_tof_pos_num = (scanner_ptr->get_max_num_timing_poss() / tof_mash_factor)-1 ;//min_tof_pos_num + (scanner_ptr->get_max_num_timing_poss() / tof_mash_factor) -1; num_tof_bins = max_tof_pos_num - min_tof_pos_num +1 ; // Ensure that we have a central tof bin. - if (num_tof_bins%2 == 0) - error("ProjDataInfo: Number of TOF bins should be an odd number. Abort."); +// if (num_tof_bins%2 == 0) +// error("ProjDataInfo: Number of TOF bins should be an odd number. Abort."); // Upper and lower boundaries of the timing poss; tof_bin_boundaries_mm.grow(min_tof_pos_num, max_tof_pos_num); diff --git a/src/buildblock/ProjDataInfoCylindrical.cxx b/src/buildblock/ProjDataInfoCylindrical.cxx index fc5ed1d7e3..3d9ee9036f 100644 --- a/src/buildblock/ProjDataInfoCylindrical.cxx +++ b/src/buildblock/ProjDataInfoCylindrical.cxx @@ -588,13 +588,13 @@ get_LOR(LORInAxialAndNoArcCorrSinogramCoordinates& lor, const float z1 = (m_in_mm - max_a*tantheta); const float z2 = (m_in_mm - min_a*tantheta); - + const bool swap = false;//get_tof_delta_time(bin)<0; lor = LORInAxialAndNoArcCorrSinogramCoordinates(z1, z2, phi, asin(s_in_mm/get_ring_radius()), get_ring_radius(), - false);// needs to set "swapped" to false given above code + swap);// needs to set "swapped" to false given above code } #if 0 diff --git a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx index 1c101e24b9..dd4e122aa9 100644 --- a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx @@ -520,19 +520,19 @@ find_cartesian_coordinates_given_scanner_coordinates (CartesianCoordinate3Dinitialise_det1det2_to_uncompressed_view_tangpos_if_not_done_yet(); if (!det1det2_to_uncompressed_view_tangpos[det1][det2].swap_detectors) - { - d1 = det2; - d2 = det1; - r1 = Ring_B; - r2 = Ring_A; - } - else { d1 = det1; d2 = det2; r1 = Ring_A; r2 = Ring_B; } + else + { + d1 = det2; + d2 = det1; + r1 = Ring_B; + r2 = Ring_A; + } #if 0 const float df1 = (2.*_PI/num_detectors_per_ring)*(det1); @@ -562,7 +562,7 @@ find_cartesian_coordinates_given_scanner_coordinates (CartesianCoordinate3D 0) { if ( do_symmetry_swap_s && s < 0) - return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq(view180, axial_pos_shift, z_shift, transform_z); + return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq(view180, axial_pos_shift, z_shift, transform_z, proj_data_info_ptr->get_scanner_ptr()->get_max_num_timing_poss()); else { if (z_shift==0) @@ -369,13 +373,13 @@ find_sym_op_general_bin( return new SymmetryOperation_PET_CartesianGrid_swap_zq(view180, axial_pos_shift, z_shift, transform_z); else*/ if ( do_symmetry_swap_s && s < 0) - return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift); + return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_scanner_ptr()->get_max_num_timing_poss()); else return new SymmetryOperation_PET_CartesianGrid_swap_zq(view180, axial_pos_shift, z_shift, transform_z); // s > 0 } else // segment_num = 0 { - if ( do_symmetry_swap_s && s < 0) return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift); + if ( do_symmetry_swap_s && s < 0) return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_scanner_ptr()->get_max_num_timing_poss()); else { if (z_shift==0) @@ -484,7 +488,7 @@ find_basic_bin(int &segment_num, int &view_num, int &axial_pos_num, int &tangent { //when swap_s, must invert timing pos for lor probs. Symmetry operation should correct bin tangential_pos_num *= -1; - timing_pos_num *= -1; + timing_pos_num = proj_data_info_ptr->get_scanner_ptr()->get_max_num_timing_poss() - timing_pos_num; change=true; } if ( do_symmetry_shift_z && axial_pos_num != 0 ) { axial_pos_num = 0; change = true; } diff --git a/src/include/stir/recon_buildblock/SymmetryOperations_PET_CartesianGrid.h b/src/include/stir/recon_buildblock/SymmetryOperations_PET_CartesianGrid.h index c61bd1fcfe..2b92a700c8 100644 --- a/src/include/stir/recon_buildblock/SymmetryOperations_PET_CartesianGrid.h +++ b/src/include/stir/recon_buildblock/SymmetryOperations_PET_CartesianGrid.h @@ -339,8 +339,8 @@ class SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq : public SymmetryOpera private: typedef SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq self; public: - SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq(const int num_views, const int axial_pos_shift, const int z_shift, const int q) - : view180(num_views), axial_pos_shift(axial_pos_shift), z_shift(z_shift), q(q) + SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq(const int num_views, const int axial_pos_shift, const int z_shift, const int q, const int max_tof_pos_num) + : view180(num_views), axial_pos_shift(axial_pos_shift), z_shift(z_shift), q(q), timing_shift(max_tof_pos_num) {} inline void @@ -364,6 +364,7 @@ class SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq : public SymmetryOpera int axial_pos_shift; int z_shift; int q; + int timing_shift; }; class SymmetryOperation_PET_CartesianGrid_swap_xy_ymx_zq : public SymmetryOperation @@ -497,8 +498,8 @@ class SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy : public SymmetryOperatio private: typedef SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy self; public: - SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(const int num_views, const int axial_pos_shift, const int z_shift) - : view180(num_views), axial_pos_shift(axial_pos_shift), z_shift(z_shift) + SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(const int num_views, const int axial_pos_shift, const int z_shift, const int max_tof_pos_num) + : view180(num_views), axial_pos_shift(axial_pos_shift), z_shift(z_shift), timing_shift(max_tof_pos_num) {} inline void @@ -521,6 +522,7 @@ class SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy : public SymmetryOperatio int view180; int axial_pos_shift; int z_shift; + int timing_shift; }; class SymmetryOperation_PET_CartesianGrid_swap_xmy_ymx_zq : public SymmetryOperation diff --git a/src/include/stir/recon_buildblock/SymmetryOperations_PET_CartesianGrid.inl b/src/include/stir/recon_buildblock/SymmetryOperations_PET_CartesianGrid.inl index a40639501b..8f339a3555 100644 --- a/src/include/stir/recon_buildblock/SymmetryOperations_PET_CartesianGrid.inl +++ b/src/include/stir/recon_buildblock/SymmetryOperations_PET_CartesianGrid.inl @@ -351,7 +351,7 @@ SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq:: { b.axial_pos_num() += axial_pos_shift; b.tangential_pos_num() *= -1; - b.timing_pos_num() *= -1; + b.timing_pos_num() = timing_shift - b.timing_pos_num(); } void @@ -553,7 +553,7 @@ SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy:: b.axial_pos_num() += axial_pos_shift; b.segment_num() *= -1; b.tangential_pos_num() *= -1; - b.timing_pos_num() *= -1; + b.timing_pos_num() = timing_shift - b.timing_pos_num(); } void SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy:: diff --git a/src/test/test_time_of_flight.cxx b/src/test/test_time_of_flight.cxx index ae53efceaf..06a7fd6723 100644 --- a/src/test/test_time_of_flight.cxx +++ b/src/test/test_time_of_flight.cxx @@ -312,7 +312,7 @@ void TOF_Tests::test_CListEventROOT() check_if_equal(det_pos_swapped.pos1(), det_pos.pos1(), "CListEventROOT: get_detection_position with swapped detectors: equal timing_pos, but different pos1"); check_if_equal(det_pos_swapped.pos1(), det_pos.pos1(), "CListEventROOT: get_detection_position with swapped detectors: equal timing_pos, but different pos1"); } - else if (det_pos_swapped.timing_pos() == -det_pos.timing_pos()) + else if (det_pos_swapped.timing_pos() == event.get_scanner_ptr()->get_max_num_timing_poss()-1-det_pos.timing_pos()) { check_if_equal(det_pos_swapped.pos2(), det_pos.pos1(), "CListEventROOT: get_detection_position with swapped detectors: equal timing_pos, but different pos1"); check_if_equal(det_pos_swapped.pos2(), det_pos.pos1(), "CListEventROOT: get_detection_position with swapped detectors: equal timing_pos, but different pos1"); From 5b9d603a623fdf670f3e6d8db292ea534725297b Mon Sep 17 00:00:00 2001 From: danieldeidda Date: Tue, 10 Oct 2023 15:55:40 +0100 Subject: [PATCH 2/9] furter adaptation to timing_pos_num from 0 to num_tof_bin -1, more swaps found --- .../ProjDataInfoCylindricalNoArcCorr.cxx | 31 ++++++++++++------- .../stir/ProjDataInfoCylindricalNoArcCorr.inl | 6 ++-- ...CListEventScannerWithDiscreteDetectors.inl | 2 +- ...ataSymmetriesForBins_PET_CartesianGrid.inl | 4 --- ...ihoodWithLinearModelForMeanAndProjData.cxx | 10 +++--- src/test/test_proj_data_info.cxx | 6 ++-- src/test/test_time_of_flight.cxx | 7 +++-- 7 files changed, 35 insertions(+), 31 deletions(-) diff --git a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx index dd4e122aa9..eafaa9f806 100644 --- a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx @@ -355,7 +355,7 @@ get_num_det_pos_pairs_for_bin(const Bin& bin) const { return get_num_ring_pairs_for_segment_axial_pos_num(bin.segment_num(), - bin.axial_pos_num())* + bin.axial_pos_num())* get_view_mashing_factor()* std::max(1,get_tof_mash_factor()); } @@ -389,24 +389,30 @@ get_all_det_pos_pairs_for_bin(vector >& dps, rings_iter != ring_pairs.end(); ++rings_iter) { - for (int uncompressed_timing_pos_num = bin.timing_pos_num()*get_tof_mash_factor() - (get_tof_mash_factor() / 2); - uncompressed_timing_pos_num <= bin.timing_pos_num()*get_tof_mash_factor() + (get_tof_mash_factor() / 2); - ++uncompressed_timing_pos_num) + int max_utpos = 0; + if (get_tof_mash_factor()==0) + max_utpos = 1; + else + max_utpos = bin.timing_pos_num()*get_tof_mash_factor() + (get_tof_mash_factor()); + + for (int uncompressed_timing_pos_num = bin.timing_pos_num()*get_tof_mash_factor() ; + uncompressed_timing_pos_num < max_utpos; + uncompressed_timing_pos_num++) { assert(current_dp_num < get_num_det_pos_pairs_for_bin(bin)); dps[current_dp_num].pos1().tangential_coord() = det1_num; dps[current_dp_num].pos1().axial_coord() = rings_iter->first; dps[current_dp_num].pos2().tangential_coord() = det2_num; dps[current_dp_num].pos2().axial_coord() = rings_iter->second; - // need to keep dp.timing_pos positive - if (uncompressed_timing_pos_num > 0) + // need to keep dp.timing_pos-get_max_tof_pos_num()/2.F positive + if (uncompressed_timing_pos_num - get_tof_mash_factor()*get_max_tof_pos_num()/2>= 0) { dps[current_dp_num].timing_pos() = static_cast(uncompressed_timing_pos_num); } else { std::swap(dps[current_dp_num].pos1(), dps[current_dp_num].pos2()); - dps[current_dp_num].timing_pos() = static_cast(-uncompressed_timing_pos_num); + dps[current_dp_num].timing_pos() = get_tof_mash_factor()*get_num_tof_poss()-static_cast(uncompressed_timing_pos_num)-1; } ++current_dp_num; } @@ -515,11 +521,11 @@ find_cartesian_coordinates_given_scanner_coordinates (CartesianCoordinate3Dget_num_detectors_per_ring(); - int d1, d2, r1, r2; + int d1, d2, r1, r2, tpos; this->initialise_det1det2_to_uncompressed_view_tangpos_if_not_done_yet(); - - if (!det1det2_to_uncompressed_view_tangpos[det1][det2].swap_detectors) + bool swap = det1det2_to_uncompressed_view_tangpos[det1][det2].swap_detectors; + if (swap) { d1 = det1; d2 = det2; @@ -532,6 +538,7 @@ find_cartesian_coordinates_given_scanner_coordinates (CartesianCoordinate3Dget_tof_mash_factor()==0 ? 0 // use timing_pos==0 in the nonTOF case - : stir::round((float)dp.timing_pos()/this->get_tof_mash_factor())); + : (int)((float)dp.timing_pos()/this->get_tof_mash_factor())); } void ProjDataInfoCylindricalNoArcCorr:: @@ -189,7 +189,7 @@ get_det_pos_pair_for_bin( //lousy work around because types don't match (short/int). TODO remove! int t1, a1, t2, a2; get_det_pair_for_bin(t1, a1, t2, a2, bin); - if (bin.timing_pos_num()>=0) + if (bin.timing_pos_num()-get_max_tof_pos_num()/2.F>=0) { dp.pos1().tangential_coord()=t1; dp.pos1().axial_coord()=a1; diff --git a/src/include/stir/listmode/CListEventScannerWithDiscreteDetectors.inl b/src/include/stir/listmode/CListEventScannerWithDiscreteDetectors.inl index 41f187c683..81ebb68e49 100644 --- a/src/include/stir/listmode/CListEventScannerWithDiscreteDetectors.inl +++ b/src/include/stir/listmode/CListEventScannerWithDiscreteDetectors.inl @@ -61,7 +61,7 @@ CListEventScannerWithDiscreteDetectors:: get_LOR() const { LORAs2Points lor; - const bool swap = this->get_delta_time() < 0.F; + const bool swap = true;//this->get_delta_time() < 0.F; // provide somewhat shorter names for the 2 coordinates, taking swap into account CartesianCoordinate3D& coord_1 = swap ? lor.p2() : lor.p1(); CartesianCoordinate3D& coord_2 = swap ? lor.p1() : lor.p2(); diff --git a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl index 89f507f2a6..67db0a532c 100644 --- a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl +++ b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl @@ -263,10 +263,6 @@ find_sym_op_general_bin( // If doing z shifts, set the axial_pos_shift to axial_pos_num, else, set to 0 const int axial_pos_shift = do_symmetry_shift_z ? axial_pos_num : 0; - // set timing shift for timing pos transformation - Bin time_bin(segment_num, view_num, axial_pos_num, 0, 0); - find_basic_bin(time_bin.segment_num(), time_bin.view_num(), time_bin.axial_pos_num(), time_bin.tangential_pos_num(), time_bin.timing_pos_num()); - const int z_shift = do_symmetry_shift_z ? num_planes_per_axial_pos[segment_num]*axial_pos_num diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx index a9ba5c88b3..eb6d2bc2c2 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx @@ -692,14 +692,14 @@ actual_compute_subset_gradient_without_penalty(TargetT& gradient, this->proj_data_sptr, subset_num, this->num_subsets, - -this->max_segment_num_to_process, + 0, this->max_segment_num_to_process, this->zero_seg0_end_planes!=0, NULL, this->additive_proj_data_sptr, this->normalisation_sptr, caching_info_ptr, - -this->max_timing_pos_num_to_process, + 0, this->max_timing_pos_num_to_process, add_sensitivity); } @@ -741,7 +741,7 @@ actual_compute_objective_function_without_penalty(const TargetT& current_estimat this->get_time_frame_definitions().get_start_time(this->get_time_frame_num()), this->get_time_frame_definitions().get_end_time(this->get_time_frame_num()), this->caching_info_ptr, - -this->max_timing_pos_num_to_process, + 0, this->max_timing_pos_num_to_process ); @@ -857,7 +857,7 @@ add_subset_sensitivity(TargetT& sensitivity, const int subset_num) const this->get_time_frame_definitions().get_start_time(this->get_time_frame_num()), this->get_time_frame_definitions().get_end_time(this->get_time_frame_num()), this->caching_info_ptr, - use_tofsens ? -this->max_timing_pos_num_to_process : 0, + use_tofsens ? 0 : 0, use_tofsens ? this->max_timing_pos_num_to_process : 0); std::transform(sensitivity.begin_all(), sensitivity.end_all(), @@ -892,7 +892,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: add_view_seg_to_sensitivity(TargetT& sensitivity, const ViewSegmentNumbers& view_seg_nums) const { - int min_timing_pos_num = use_tofsens ? -this->max_timing_pos_num_to_process : 0; + int min_timing_pos_num = use_tofsens ? 0 : 0; int max_timing_pos_num = use_tofsens ? this->max_timing_pos_num_to_process : 0; for (int timing_pos_num = min_timing_pos_num; timing_pos_num <= max_timing_pos_num; ++timing_pos_num) { diff --git a/src/test/test_proj_data_info.cxx b/src/test/test_proj_data_info.cxx index 2d6ffdd564..287bab1a04 100644 --- a/src/test/test_proj_data_info.cxx +++ b/src/test/test_proj_data_info.cxx @@ -147,7 +147,7 @@ ProjDataInfoTests::test_generic_proj_data_info(ProjDataInfo& proj_data_info) "check on min/max_tangential_pos_num being (almost) centred"); check_if_equal(proj_data_info.get_max_tof_pos_num() - proj_data_info.get_min_tof_pos_num() + 1, proj_data_info.get_num_tof_poss(), "basic check on min/max/num_tof_pos_num"); - check_if_equal(proj_data_info.get_max_tof_pos_num() + proj_data_info.get_min_tof_pos_num(), 0, + check_if_equal(proj_data_info.get_max_tof_pos_num() + proj_data_info.get_min_tof_pos_num(), proj_data_info.get_max_tof_pos_num(), "check on min/max_tof_pos_num being (almost) centred"); check_if_equal(proj_data_info.get_max_view_num() - proj_data_info.get_min_view_num() + 1, proj_data_info.get_num_views(), "basic check on min/max/num_view_num"); @@ -236,7 +236,7 @@ ProjDataInfoTests::test_generic_proj_data_info(ProjDataInfo& proj_data_info) || !check(diff_view_num <= 1, "round-trip get_LOR then get_bin: view") || !check(diff_axial_pos_num <= 1, "round-trip get_LOR then get_bin: axial_pos") || !check(diff_tangential_pos_num <= 1, "round-trip get_LOR then get_bin: tangential_pos") - || !check(diff_timing_pos_num == 0, "round-trip get_LOR then get_bin: timing_pos")) + || !check(diff_timing_pos_num == 0 || diff_timing_pos_num == proj_data_info.get_max_tof_pos_num() , "round-trip get_LOR then get_bin: timing_pos")) #else if (!check(org_bin == new_bin, "round-trip get_LOR then get_bin")) @@ -293,7 +293,7 @@ ProjDataInfoTests::test_generic_proj_data_info(ProjDataInfo& proj_data_info) || !check(diff_view_num <= 1, "round-trip get_LOR then get_bin (LORAs2Points): view") || !check(diff_axial_pos_num <= 1, "round-trip get_LOR then get_bin (LORAs2Points): axial_pos") || !check(diff_tangential_pos_num <= 1, "round-trip get_LOR then get_bin (LORAs2Points): tangential_pos") - || !check(diff_timing_pos_num == 0, "round-trip get_LOR then get_bin: timing_pos")) + || !check(diff_timing_pos_num == 0 || diff_timing_pos_num == proj_data_info.get_max_tof_pos_num(), "round-trip get_LOR then get_bin: timing_pos")) #else if (!check(org_bin == new_bin, "round-trip get_LOR then get_bin")) diff --git a/src/test/test_time_of_flight.cxx b/src/test/test_time_of_flight.cxx index 06a7fd6723..61a3dc2b16 100644 --- a/src/test/test_time_of_flight.cxx +++ b/src/test/test_time_of_flight.cxx @@ -291,6 +291,7 @@ void TOF_Tests::test_CListEventROOT() DetectionPositionPair<> det_pos; event.get_detection_position(det_pos); + LORAs2Points lor_2pts(event.get_LOR()); LORInAxialAndNoArcCorrSinogramCoordinates lor_sc; test_proj_data_info_sptr->get_LOR(lor_sc, bin); @@ -310,12 +311,12 @@ void TOF_Tests::test_CListEventROOT() if (det_pos_swapped.timing_pos() == det_pos.timing_pos()) { check_if_equal(det_pos_swapped.pos1(), det_pos.pos1(), "CListEventROOT: get_detection_position with swapped detectors: equal timing_pos, but different pos1"); - check_if_equal(det_pos_swapped.pos1(), det_pos.pos1(), "CListEventROOT: get_detection_position with swapped detectors: equal timing_pos, but different pos1"); + check_if_equal(det_pos_swapped.pos2(), det_pos.pos2(), "CListEventROOT: get_detection_position with swapped detectors: equal timing_pos, but different pos2"); } else if (det_pos_swapped.timing_pos() == event.get_scanner_ptr()->get_max_num_timing_poss()-1-det_pos.timing_pos()) { - check_if_equal(det_pos_swapped.pos2(), det_pos.pos1(), "CListEventROOT: get_detection_position with swapped detectors: equal timing_pos, but different pos1"); - check_if_equal(det_pos_swapped.pos2(), det_pos.pos1(), "CListEventROOT: get_detection_position with swapped detectors: equal timing_pos, but different pos1"); + check_if_equal(det_pos_swapped.pos2(), det_pos.pos1(), "CListEventROOT: get_detection_position with swapped detectors: opposite timing_pos, but different pos1/2"); + check_if_equal(det_pos_swapped.pos1(), det_pos.pos2(), "CListEventROOT: get_detection_position with swapped detectors: opposite timing_pos, but different pos1/2"); } else { From da07d940a81d50a38ee7aee61df73ce40dff70b9 Mon Sep 17 00:00:00 2001 From: danieldeidda Date: Tue, 17 Oct 2023 14:33:57 +0100 Subject: [PATCH 3/9] furter adaptation to timing_pos_num from 0 to num_tof_bin -1, fixed get_k and swap s sym --- src/buildblock/ProjDataInfo.cxx | 5 +---- .../DataSymmetriesForBins_PET_CartesianGrid.inl | 2 +- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/src/buildblock/ProjDataInfo.cxx b/src/buildblock/ProjDataInfo.cxx index cf391c0447..f78f14e17d 100644 --- a/src/buildblock/ProjDataInfo.cxx +++ b/src/buildblock/ProjDataInfo.cxx @@ -73,10 +73,7 @@ START_NAMESPACE_STIR float ProjDataInfo::get_k(const Bin& bin) const { - if (!(num_tof_bins%2)) - return (bin.timing_pos_num() - (num_tof_bins)/2.f)* tof_increament_in_mm ; - else - return (bin.timing_pos_num() - (num_tof_bins-1)/2.f)* tof_increament_in_mm ; + return (bin.timing_pos_num() - (num_tof_bins-1)/2.f)* tof_increament_in_mm ; } double diff --git a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl index 67db0a532c..e1d1ebf1f1 100644 --- a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl +++ b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl @@ -484,7 +484,7 @@ find_basic_bin(int &segment_num, int &view_num, int &axial_pos_num, int &tangent { //when swap_s, must invert timing pos for lor probs. Symmetry operation should correct bin tangential_pos_num *= -1; - timing_pos_num = proj_data_info_ptr->get_scanner_ptr()->get_max_num_timing_poss() - timing_pos_num; + timing_pos_num = proj_data_info_ptr->get_max_tof_pos_num() - timing_pos_num; change=true; } if ( do_symmetry_shift_z && axial_pos_num != 0 ) { axial_pos_num = 0; change = true; } From e9ffdd43221ca5737d1125bf9b4dcd4222febb87 Mon Sep 17 00:00:00 2001 From: danieldeidda Date: Mon, 15 Jan 2024 15:54:47 +0000 Subject: [PATCH 4/9] more tof swaps and fixes for even tof bins --- src/buildblock/ProjDataInfo.cxx | 2 +- .../ProjDataInfoCylindricalNoArcCorr.cxx | 34 ++++++++---------- ...elihoodWithLinearModelForMeanAndProjData.h | 7 ++-- ...ihoodWithLinearModelForMeanAndProjData.cxx | 36 ++++++++++--------- 4 files changed, 39 insertions(+), 40 deletions(-) diff --git a/src/buildblock/ProjDataInfo.cxx b/src/buildblock/ProjDataInfo.cxx index f78f14e17d..4191b41a40 100644 --- a/src/buildblock/ProjDataInfo.cxx +++ b/src/buildblock/ProjDataInfo.cxx @@ -73,7 +73,7 @@ START_NAMESPACE_STIR float ProjDataInfo::get_k(const Bin& bin) const { - return (bin.timing_pos_num() - (num_tof_bins-1)/2.f)* tof_increament_in_mm ; + return (bin.timing_pos_num() - (num_tof_bins-1)/2.f)* tof_increament_in_mm ; } double diff --git a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx index eafaa9f806..50bec2771b 100644 --- a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx @@ -404,16 +404,9 @@ get_all_det_pos_pairs_for_bin(vector >& dps, dps[current_dp_num].pos1().axial_coord() = rings_iter->first; dps[current_dp_num].pos2().tangential_coord() = det2_num; dps[current_dp_num].pos2().axial_coord() = rings_iter->second; - // need to keep dp.timing_pos-get_max_tof_pos_num()/2.F positive - if (uncompressed_timing_pos_num - get_tof_mash_factor()*get_max_tof_pos_num()/2>= 0) - { - dps[current_dp_num].timing_pos() = static_cast(uncompressed_timing_pos_num); - } - else - { - std::swap(dps[current_dp_num].pos1(), dps[current_dp_num].pos2()); - dps[current_dp_num].timing_pos() = get_tof_mash_factor()*get_num_tof_poss()-static_cast(uncompressed_timing_pos_num)-1; - } + + dps[current_dp_num].timing_pos() = static_cast(uncompressed_timing_pos_num); + ++current_dp_num; } } @@ -524,15 +517,8 @@ find_cartesian_coordinates_given_scanner_coordinates (CartesianCoordinate3Dinitialise_det1det2_to_uncompressed_view_tangpos_if_not_done_yet(); - bool swap = det1det2_to_uncompressed_view_tangpos[det1][det2].swap_detectors; - if (swap) - { - d1 = det1; - d2 = det2; - r1 = Ring_A; - r2 = Ring_B; - } - else + + if (!det1det2_to_uncompressed_view_tangpos[det1][det2].swap_detectors) { d1 = det2; d2 = det1; @@ -540,6 +526,14 @@ find_cartesian_coordinates_given_scanner_coordinates (CartesianCoordinate3D& get_proj_data_sptr() const; const int get_max_segment_num_to_process() const; - const int get_max_timing_pos_num_to_process() const; + const int get_max_num_central_timing_poss_to_process() const; const bool get_zero_seg0_end_planes() const; const ProjData& get_additive_proj_data() const; const shared_ptr& get_additive_proj_data_sptr() const; @@ -201,7 +201,7 @@ public RegisteredParsingObject&); void set_max_segment_num_to_process(const int); - void set_max_timing_pos_num_to_process(const int); + void set_max_num_central_timing_poss_to_process(const int); void set_zero_seg0_end_planes(const bool); //N.E. Changed to ExamData virtual void set_additive_proj_data_sptr(const shared_ptr&); @@ -311,7 +311,8 @@ public RegisteredParsingObject target_parameter_parser; diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx index eb6d2bc2c2..a5aa34ed3a 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx @@ -97,7 +97,7 @@ set_defaults() this->input_filename=""; this->max_segment_num_to_process=-1; - this->max_timing_pos_num_to_process=0; + this->max_num_central_timing_poss_to_process=0; // KT 20/06/2001 disabled //num_views_to_add=1; this->proj_data_sptr.reset(); //MJ added @@ -324,8 +324,8 @@ get_max_segment_num_to_process() const template const int PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: -get_max_timing_pos_num_to_process() const -{ return this->max_timing_pos_num_to_process; } +get_max_num_central_timing_poss_to_process() const +{ return this->max_num_central_timing_poss_to_process; } template const bool @@ -415,9 +415,9 @@ set_max_segment_num_to_process(const int arg) template void PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: -set_max_timing_pos_num_to_process(const int arg) +set_max_num_central_timing_poss_to_process(const int arg) { - this->max_timing_pos_num_to_process = arg; + this->max_num_central_timing_poss_to_process = arg; } template @@ -569,8 +569,11 @@ set_up_before_sensitivity(shared_ptr const& target_sptr) return Succeeded::no; } - this->max_timing_pos_num_to_process = + this->max_num_central_timing_poss_to_process = this->proj_data_sptr->get_max_tof_pos_num(); +// The following "shift"is used to estimate the maximum and minimum tof position number to process +// i.e. min_pos= num_tof_bins+shift and max_pos= num_tof_bins-shift if we use all than shift =0 + this->tof_bin_shift = (proj_data_sptr->get_num_tof_poss()-this->max_num_central_timing_poss_to_process)/2; shared_ptr proj_data_info_sptr(this->proj_data_sptr->get_proj_data_info_sptr()->clone()); @@ -699,8 +702,8 @@ actual_compute_subset_gradient_without_penalty(TargetT& gradient, this->additive_proj_data_sptr, this->normalisation_sptr, caching_info_ptr, - 0, - this->max_timing_pos_num_to_process, + proj_data_sptr->get_min_tof_pos_num() + tof_bin_shift, + proj_data_sptr->get_max_tof_pos_num() - tof_bin_shift, add_sensitivity); } @@ -741,8 +744,8 @@ actual_compute_objective_function_without_penalty(const TargetT& current_estimat this->get_time_frame_definitions().get_start_time(this->get_time_frame_num()), this->get_time_frame_definitions().get_end_time(this->get_time_frame_num()), this->caching_info_ptr, - 0, - this->max_timing_pos_num_to_process + proj_data_sptr->get_min_tof_pos_num() + tof_bin_shift, + proj_data_sptr->get_max_tof_pos_num() - tof_bin_shift ); @@ -757,10 +760,11 @@ sum_projection_data() const { float counts=0.0F; - + int min_timing_poss_to_process = proj_data_sptr->get_min_tof_pos_num() + tof_bin_shift; + int max_timing_poss_to_process = proj_data_sptr->get_max_tof_pos_num() - tof_bin_shift; for (int segment_num = -max_segment_num_to_process; segment_num <= max_segment_num_to_process; ++segment_num) { - for (int timing_pos_num = -max_timing_pos_num_to_process; timing_pos_num <= max_timing_pos_num_to_process; ++timing_pos_num) + for (int timing_pos_num = min_timing_poss_to_process; timing_pos_num <= max_timing_poss_to_process; ++timing_pos_num) { for (int view_num = proj_data_sptr->get_min_view_num(); view_num <= proj_data_sptr->get_max_view_num(); @@ -857,8 +861,8 @@ add_subset_sensitivity(TargetT& sensitivity, const int subset_num) const this->get_time_frame_definitions().get_start_time(this->get_time_frame_num()), this->get_time_frame_definitions().get_end_time(this->get_time_frame_num()), this->caching_info_ptr, - use_tofsens ? 0 : 0, - use_tofsens ? this->max_timing_pos_num_to_process : 0); + use_tofsens ? proj_data_sptr->get_min_tof_pos_num() + tof_bin_shift : 0, + use_tofsens ? proj_data_sptr->get_max_tof_pos_num() - tof_bin_shift : 0); std::transform(sensitivity.begin_all(), sensitivity.end_all(), sensitivity_this_subset_sptr->begin_all(), sensitivity.begin_all(), @@ -892,8 +896,8 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: add_view_seg_to_sensitivity(TargetT& sensitivity, const ViewSegmentNumbers& view_seg_nums) const { - int min_timing_pos_num = use_tofsens ? 0 : 0; - int max_timing_pos_num = use_tofsens ? this->max_timing_pos_num_to_process : 0; + int min_timing_pos_num = use_tofsens ? proj_data_sptr->get_min_tof_pos_num() + tof_bin_shift : 0; + int max_timing_pos_num = use_tofsens ? proj_data_sptr->get_max_tof_pos_num() - tof_bin_shift : 0; for (int timing_pos_num = min_timing_pos_num; timing_pos_num <= max_timing_pos_num; ++timing_pos_num) { RelatedViewgrams viewgrams = From 605f950f56b6536a3dbe16658d9389219ebec1f2 Mon Sep 17 00:00:00 2001 From: danieldeidda Date: Mon, 22 Jan 2024 16:54:49 +0000 Subject: [PATCH 5/9] fixing symmetry test issues --- .../DataSymmetriesForBins_PET_CartesianGrid.inl | 6 +++--- .../test_DataSymmetriesForBins_PET_CartesianGrid.cxx | 12 ++---------- 2 files changed, 5 insertions(+), 13 deletions(-) diff --git a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl index e1d1ebf1f1..1bea91c3d1 100644 --- a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl +++ b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl @@ -353,7 +353,7 @@ find_sym_op_general_bin( if ( !do_symmetry_swap_segment || segment_num > 0) { if ( do_symmetry_swap_s && s < 0) - return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq(view180, axial_pos_shift, z_shift, transform_z, proj_data_info_ptr->get_scanner_ptr()->get_max_num_timing_poss()); + return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq(view180, axial_pos_shift, z_shift, transform_z, proj_data_info_ptr->get_scanner_ptr()->get_max_num_timing_poss()-1); else { if (z_shift==0) @@ -369,13 +369,13 @@ find_sym_op_general_bin( return new SymmetryOperation_PET_CartesianGrid_swap_zq(view180, axial_pos_shift, z_shift, transform_z); else*/ if ( do_symmetry_swap_s && s < 0) - return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_scanner_ptr()->get_max_num_timing_poss()); + return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_scanner_ptr()->get_max_num_timing_poss()-1); else return new SymmetryOperation_PET_CartesianGrid_swap_zq(view180, axial_pos_shift, z_shift, transform_z); // s > 0 } else // segment_num = 0 { - if ( do_symmetry_swap_s && s < 0) return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_scanner_ptr()->get_max_num_timing_poss()); + if ( do_symmetry_swap_s && s < 0) return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_scanner_ptr()->get_max_num_timing_poss()-1); else { if (z_shift==0) diff --git a/src/recon_test/test_DataSymmetriesForBins_PET_CartesianGrid.cxx b/src/recon_test/test_DataSymmetriesForBins_PET_CartesianGrid.cxx index 24dd7790b0..1195bc799a 100644 --- a/src/recon_test/test_DataSymmetriesForBins_PET_CartesianGrid.cxx +++ b/src/recon_test/test_DataSymmetriesForBins_PET_CartesianGrid.cxx @@ -130,16 +130,8 @@ run_tests_2_proj_matrices_1_bin(const ProjMatrixByBin& proj_matrix_no_symm, { // SYM const Bin bin=*bin_iter; - cerr << "\nCurrent bin: \tsegment = " << bin.segment_num() - << ", \taxial pos " << bin.axial_pos_num() - << ", \tview = " << bin.view_num() - << ", \ttangential_pos_num = " << bin.tangential_pos_num() - << ", timing position index = " << bin.timing_pos_num() - << "\nSymm bin: \t\tsegment = " << elems_with_sym.get_bin().segment_num() - << ", \taxial pos " << elems_with_sym.get_bin().axial_pos_num() - << ", \tview = " << elems_with_sym.get_bin().view_num() - << ", \ttangential_pos_num = " << elems_with_sym.get_bin().tangential_pos_num() - << ", timing position index = " << bin.timing_pos_num() << "\n"; + cerr << "\nCurrent bin: " << elems_no_sym.get_bin() + << "\nSymm bin: " << elems_with_sym.get_bin()<< "\n"; if (elems_no_sym != elems_with_sym) { From fcfb1caec4e741cc1b9a31ce0e4d691598dd7332 Mon Sep 17 00:00:00 2001 From: danieldeidda Date: Fri, 26 Jan 2024 15:35:02 +0000 Subject: [PATCH 6/9] fixing asserts and tests --- .../ProjDataInfoCylindricalNoArcCorr.cxx | 5 +++- .../stir/ProjDataInfoCylindricalNoArcCorr.inl | 2 +- ...ihoodWithLinearModelForMeanAndProjData.cxx | 4 ++-- src/test/test_proj_data_info.cxx | 24 ++++++++++--------- 4 files changed, 20 insertions(+), 15 deletions(-) diff --git a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx index f503d45d04..c1d3d191a8 100644 --- a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx @@ -369,7 +369,10 @@ get_all_det_pos_pairs_for_bin(vector >& dps, assert(!is_tof_data() || (get_tof_mash_factor() % 2 == 1)); // TODOTOF // we will need to add all (unmashed) timing_pos for the current bin min_timing_pos_num = bin.timing_pos_num()*get_tof_mash_factor();//bin.timing_pos_num()*get_tof_mash_factor() - (get_tof_mash_factor() / 2); - max_timing_pos_num = bin.timing_pos_num()*get_tof_mash_factor() + (get_tof_mash_factor());//bin.timing_pos_num()*get_tof_mash_factor() + (get_tof_mash_factor() / 2); + if (get_tof_mash_factor()==0) + max_timing_pos_num = bin.timing_pos_num();//bin.timing_pos_num()*get_tof_mash_factor() + (get_tof_mash_factor() / 2); + else + max_timing_pos_num = bin.timing_pos_num()*get_tof_mash_factor() + (get_tof_mash_factor())-1; } diff --git a/src/include/stir/ProjDataInfoCylindricalNoArcCorr.inl b/src/include/stir/ProjDataInfoCylindricalNoArcCorr.inl index b645e5f492..7acb78c4bc 100644 --- a/src/include/stir/ProjDataInfoCylindricalNoArcCorr.inl +++ b/src/include/stir/ProjDataInfoCylindricalNoArcCorr.inl @@ -156,7 +156,7 @@ get_bin_for_det_pos_pair(Bin& bin, dp.pos2().axial_coord(), this->get_tof_mash_factor()==0 ? 0 // use timing_pos==0 in the nonTOF case - : (int)((float)dp.timing_pos()/this->get_tof_mash_factor())); + : (int)(dp.timing_pos()/this->get_tof_mash_factor())); } void ProjDataInfoCylindricalNoArcCorr:: diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx index e0981e4dad..a6a07be90e 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx @@ -97,7 +97,7 @@ set_defaults() this->input_filename=""; this->max_segment_num_to_process=-1; - this->max_num_central_timing_poss_to_process=0; + this->max_num_central_timing_poss_to_process=-1; // KT 20/06/2001 disabled //num_views_to_add=1; this->proj_data_sptr.reset(); //MJ added @@ -732,7 +732,7 @@ actual_compute_subset_gradient_without_penalty(TargetT& gradient, this->proj_data_sptr, subset_num, this->num_subsets, - 0, + -this->max_segment_num_to_process, this->max_segment_num_to_process, this->zero_seg0_end_planes!=0, NULL, diff --git a/src/test/test_proj_data_info.cxx b/src/test/test_proj_data_info.cxx index 436ebd9ff5..8abea1c9da 100644 --- a/src/test/test_proj_data_info.cxx +++ b/src/test/test_proj_data_info.cxx @@ -147,7 +147,7 @@ ProjDataInfoTests::test_generic_proj_data_info(ProjDataInfo& proj_data_info) "check on min/max_tangential_pos_num being (almost) centred"); check_if_equal(proj_data_info.get_max_tof_pos_num() - proj_data_info.get_min_tof_pos_num() + 1, proj_data_info.get_num_tof_poss(), "basic check on min/max/num_tof_pos_num"); - check_if_equal(proj_data_info.get_max_tof_pos_num() + proj_data_info.get_min_tof_pos_num(), proj_data_info.get_max_tof_pos_num(), + check_if_equal(proj_data_info.get_min_tof_pos_num() , 0, "check on min/max_tof_pos_num being (almost) centred"); check_if_equal(proj_data_info.get_max_view_num() - proj_data_info.get_min_view_num() + 1, proj_data_info.get_num_views(), "basic check on min/max/num_view_num"); @@ -204,16 +204,17 @@ ProjDataInfoTests::test_generic_proj_data_info(ProjDataInfo& proj_data_info) { const Bin new_bin = proj_data_info.get_bin(lor, delta_time); + const bool views_are_close = intabs(org_bin.view_num() - new_bin.view_num()) < proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num()); #if 1 // the differences need to also consider wrap-around in views, which would flip tangential pos and segment and TOF bin - const int diff_segment_num = intabs(org_bin.view_num() - new_bin.view_num()) < proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num()) ? + const int diff_segment_num = views_are_close ? intabs(org_bin.segment_num() - new_bin.segment_num()) : intabs(org_bin.segment_num() + new_bin.segment_num()); const int diff_view_num = min(intabs(org_bin.view_num() - new_bin.view_num()), proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num())); const int diff_axial_pos_num = intabs(org_bin.axial_pos_num() - new_bin.axial_pos_num()); - const int diff_tangential_pos_num = intabs(org_bin.view_num() - new_bin.view_num()) < proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num()) ? + const int diff_tangential_pos_num = views_are_close ? intabs(org_bin.tangential_pos_num() - new_bin.tangential_pos_num()) : intabs(org_bin.tangential_pos_num() + new_bin.tangential_pos_num()); - const int diff_timing_pos_num = intabs(org_bin.view_num() - new_bin.view_num()) < proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num()) ? - intabs(org_bin.timing_pos_num() - new_bin.timing_pos_num()) : intabs(org_bin.timing_pos_num() + new_bin.timing_pos_num()); + const int diff_timing_pos_num = views_are_close ? + intabs(org_bin.timing_pos_num() - new_bin.timing_pos_num()) : intabs(org_bin.timing_pos_num() - (proj_data_info.get_max_tof_pos_num() - new_bin.timing_pos_num())); if (new_bin.get_bin_value() > 0) { if (diff_segment_num > max_diff_segment_num) @@ -236,7 +237,7 @@ ProjDataInfoTests::test_generic_proj_data_info(ProjDataInfo& proj_data_info) || !check(diff_view_num <= 1, "round-trip get_LOR then get_bin: view") || !check(diff_axial_pos_num <= 1, "round-trip get_LOR then get_bin: axial_pos") || !check(diff_tangential_pos_num <= 1, "round-trip get_LOR then get_bin: tangential_pos") - || !check(diff_timing_pos_num == 0 || diff_timing_pos_num == proj_data_info.get_max_tof_pos_num() , "round-trip get_LOR then get_bin: timing_pos")) + || !check(diff_timing_pos_num == 0 , "round-trip get_LOR then get_bin: timing_pos")) #else if (!check(org_bin == new_bin, "round-trip get_LOR then get_bin")) @@ -263,17 +264,18 @@ ProjDataInfoTests::test_generic_proj_data_info(ProjDataInfo& proj_data_info) << ", y2=" << lor_as_points.p2().y() << ", x2=" << lor_as_points.p2().x() << std::endl; #endif const Bin new_bin = proj_data_info.get_bin(lor_as_points, proj_data_info.get_tof_delta_time(org_bin)); + const bool views_are_close = intabs(org_bin.view_num() - new_bin.view_num()) < proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num()); #if 1 // the differences need to also consider wrap-around in views, which would flip tangential pos and segment - const int diff_segment_num = intabs(org_bin.view_num() - new_bin.view_num()) < proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num()) ? + const int diff_segment_num = views_are_close ? intabs(org_bin.segment_num() - new_bin.segment_num()) : intabs(org_bin.segment_num() + new_bin.segment_num()); const int diff_view_num = min(intabs(org_bin.view_num() - new_bin.view_num()), proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num())); const int diff_axial_pos_num = intabs(org_bin.axial_pos_num() - new_bin.axial_pos_num()); - const int diff_tangential_pos_num = intabs(org_bin.view_num() - new_bin.view_num()) < proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num()) ? + const int diff_tangential_pos_num = views_are_close ? intabs(org_bin.tangential_pos_num() - new_bin.tangential_pos_num()) : intabs(org_bin.tangential_pos_num() + new_bin.tangential_pos_num()); - const int diff_timing_pos_num = intabs(org_bin.view_num() - new_bin.view_num()) < proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num()) ? - intabs(org_bin.timing_pos_num() - new_bin.timing_pos_num()) : intabs(org_bin.timing_pos_num() + new_bin.timing_pos_num()); + const int diff_timing_pos_num = views_are_close ? + intabs(org_bin.timing_pos_num() - new_bin.timing_pos_num()) : intabs(org_bin.timing_pos_num() - (proj_data_info.get_max_tof_pos_num() - new_bin.timing_pos_num())); if (new_bin.get_bin_value() > 0) { if (diff_segment_num > max_diff_segment_num) @@ -293,7 +295,7 @@ ProjDataInfoTests::test_generic_proj_data_info(ProjDataInfo& proj_data_info) || !check(diff_view_num <= 1, "round-trip get_LOR then get_bin (LORAs2Points): view") || !check(diff_axial_pos_num <= 1, "round-trip get_LOR then get_bin (LORAs2Points): axial_pos") || !check(diff_tangential_pos_num <= 1, "round-trip get_LOR then get_bin (LORAs2Points): tangential_pos") - || !check(diff_timing_pos_num == 0 || diff_timing_pos_num == proj_data_info.get_max_tof_pos_num(), "round-trip get_LOR then get_bin: timing_pos")) + || !check(diff_timing_pos_num == 0, "round-trip get_LOR then get_bin: timing_pos")) #else if (!check(org_bin == new_bin, "round-trip get_LOR then get_bin")) From ab7c5df8f02eaabe28193a8ce8459e3898cc6a0a Mon Sep 17 00:00:00 2001 From: danieldeidda Date: Mon, 11 Mar 2024 13:31:55 +0000 Subject: [PATCH 7/9] fixing symmetry indeces --- .../DataSymmetriesForBins_PET_CartesianGrid.inl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl index 1bea91c3d1..bfe302882a 100644 --- a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl +++ b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl @@ -353,7 +353,7 @@ find_sym_op_general_bin( if ( !do_symmetry_swap_segment || segment_num > 0) { if ( do_symmetry_swap_s && s < 0) - return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq(view180, axial_pos_shift, z_shift, transform_z, proj_data_info_ptr->get_scanner_ptr()->get_max_num_timing_poss()-1); + return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq(view180, axial_pos_shift, z_shift, transform_z, proj_data_info_ptr->get_max_tof_pos_num()-1); else { if (z_shift==0) @@ -369,13 +369,13 @@ find_sym_op_general_bin( return new SymmetryOperation_PET_CartesianGrid_swap_zq(view180, axial_pos_shift, z_shift, transform_z); else*/ if ( do_symmetry_swap_s && s < 0) - return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_scanner_ptr()->get_max_num_timing_poss()-1); + return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_max_tof_pos_num()-1); else return new SymmetryOperation_PET_CartesianGrid_swap_zq(view180, axial_pos_shift, z_shift, transform_z); // s > 0 } else // segment_num = 0 { - if ( do_symmetry_swap_s && s < 0) return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_scanner_ptr()->get_max_num_timing_poss()-1); + if ( do_symmetry_swap_s && s < 0) return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_max_tof_pos_num()-1); else { if (z_shift==0) From ab23474bd08be5ff36c712adee489a4f5a3688dc Mon Sep 17 00:00:00 2001 From: danieldeidda Date: Mon, 25 Mar 2024 17:23:24 +0000 Subject: [PATCH 8/9] fix Datasymmetry timeshift; change IO timing pos; add check in ProjDatafromStream; add glossary for TOF number --- documentation/STIR-glossary.tex | 3 +++ src/IO/stir_ecat_common.cxx | 2 +- src/buildblock/ProjDataFromStream.cxx | 4 ++++ .../DataSymmetriesForBins_PET_CartesianGrid.inl | 6 +++--- 4 files changed, 11 insertions(+), 4 deletions(-) diff --git a/documentation/STIR-glossary.tex b/documentation/STIR-glossary.tex index b4786d0d79..0aa5e7ad37 100644 --- a/documentation/STIR-glossary.tex +++ b/documentation/STIR-glossary.tex @@ -135,6 +135,9 @@ \section*{Basics} Set of \textbf{merged} \textbf{sinograms} with a common average \textbf{ring difference} as shown in Annex 1. +\item[TOF number] +The line of response is divided $N_t$ TOF positions the numbering goes from 0 to $N_t -1$ + \item[Viewgram] Set of equal azimuth \textbf{merged} \textbf{LORs} of a \textbf{segment}. diff --git a/src/IO/stir_ecat_common.cxx b/src/IO/stir_ecat_common.cxx index 0fa9df608c..e7ce585cf3 100644 --- a/src/IO/stir_ecat_common.cxx +++ b/src/IO/stir_ecat_common.cxx @@ -208,7 +208,7 @@ find_timing_poss_sequence(const ProjDataInfo& pdi) for (int timing_pos_num = 1; timing_pos_num <= max_timing_pos_num; ++timing_pos_num) { timing_pos_sequence[2 * timing_pos_num - 1] = timing_pos_num; - timing_pos_sequence[2 * timing_pos_num] = -timing_pos_num; + timing_pos_sequence[2 * timing_pos_num] = max_timing_pos_num - timing_pos_num; } return timing_pos_sequence; } diff --git a/src/buildblock/ProjDataFromStream.cxx b/src/buildblock/ProjDataFromStream.cxx index bcb47f9c65..8411f09dbd 100644 --- a/src/buildblock/ProjDataFromStream.cxx +++ b/src/buildblock/ProjDataFromStream.cxx @@ -151,6 +151,10 @@ ProjDataFromStream::activate_TOF() void ProjDataFromStream::set_timing_poss_sequence_in_stream(const std::vector& seq) { + if(!(seq[0]==this->proj_data_info_sptr->get_min_tof_pos_num() && + seq[seq.size()-1]==this->proj_data_info_sptr->get_max_tof_pos_num())) + error("the timing position sequence is different from STIR convention 0-max!"); + this->timing_poss_sequence = seq; } diff --git a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl index 85077daaca..de459093c1 100644 --- a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl +++ b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl @@ -327,7 +327,7 @@ DataSymmetriesForBins_PET_CartesianGrid::find_sym_op_general_bin(int s, int segm { if (do_symmetry_swap_s && s < 0) return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq( - view180, axial_pos_shift, z_shift, transform_z, proj_data_info_ptr->get_max_tof_pos_num() - 1); + view180, axial_pos_shift, z_shift, transform_z, proj_data_info_ptr->get_max_tof_pos_num()); else { if (z_shift == 0) @@ -343,7 +343,7 @@ DataSymmetriesForBins_PET_CartesianGrid::find_sym_op_general_bin(int s, int segm else*/ if (do_symmetry_swap_s && s < 0) return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy( - view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_max_tof_pos_num() - 1); + view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_max_tof_pos_num()); else return new SymmetryOperation_PET_CartesianGrid_swap_zq(view180, axial_pos_shift, z_shift, transform_z); // s > 0 } @@ -352,7 +352,7 @@ DataSymmetriesForBins_PET_CartesianGrid::find_sym_op_general_bin(int s, int segm if (do_symmetry_swap_s && s < 0) return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy( - view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_max_tof_pos_num() - 1); + view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_max_tof_pos_num()); else { From ab7265dd60b6d74e1ca571975b5f8078e31ca91d Mon Sep 17 00:00:00 2001 From: Daniel Deidda Date: Tue, 26 Mar 2024 09:32:46 +0000 Subject: [PATCH 9/9] Update src/IO/stir_ecat_common.cxx Co-authored-by: Kris Thielemans --- src/IO/stir_ecat_common.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/IO/stir_ecat_common.cxx b/src/IO/stir_ecat_common.cxx index e7ce585cf3..ebd49831e9 100644 --- a/src/IO/stir_ecat_common.cxx +++ b/src/IO/stir_ecat_common.cxx @@ -207,7 +207,7 @@ find_timing_poss_sequence(const ProjDataInfo& pdi) timing_pos_sequence[0] = 0; for (int timing_pos_num = 1; timing_pos_num <= max_timing_pos_num; ++timing_pos_num) { - timing_pos_sequence[2 * timing_pos_num - 1] = timing_pos_num; + timing_pos_sequence[2 * timing_pos_num - 1] = max_timing_pos_num + timing_pos_num; timing_pos_sequence[2 * timing_pos_num] = max_timing_pos_num - timing_pos_num; } return timing_pos_sequence;