diff --git a/src/buildblock/extend_projdata.cxx b/src/buildblock/extend_projdata.cxx index ae9496f01e..33df31bcab 100644 --- a/src/buildblock/extend_projdata.cxx +++ b/src/buildblock/extend_projdata.cxx @@ -119,10 +119,11 @@ namespace detail Array<3,float> extend_segment_in_views(const SegmentBySinogram& sino, - const int min_view_extension, const int max_view_extension) + const int min_view_extension, const int max_view_extension, + const SegmentBySinogram& inv_sino) { - if (sino.get_segment_num()!=0) - error("extend_segment with single segment works only for segment 0"); +// if (sino.get_segment_num()!=0) +// error("extend_segment with single segment works only for segment 0"); BasicCoordinate<3,int> min, max; @@ -138,7 +139,7 @@ extend_segment_in_views(const SegmentBySinogram& sino, { out[ax_pos_num] = detail:: - extend_sinogram_in_views(sino[ax_pos_num],sino[ax_pos_num], + extend_sinogram_in_views(sino[ax_pos_num],inv_sino[ax_pos_num], *(sino.get_proj_data_info_sptr()), min_view_extension, max_view_extension); } @@ -147,14 +148,15 @@ extend_segment_in_views(const SegmentBySinogram& sino, Array<2,float> extend_sinogram_in_views(const Sinogram& sino, - const int min_view_extension, const int max_view_extension) + const int min_view_extension, const int max_view_extension, + const Sinogram& inv_sino) { - if (sino.get_segment_num()!=0) - error("extend_segment with single segment works only for segment 0"); +// if (sino.get_segment_num()!=0) +// error("extend_segment with single segment works only for segment 0"); return detail:: - extend_sinogram_in_views(sino, sino, + extend_sinogram_in_views(sino, inv_sino, *(sino.get_proj_data_info_sptr()), min_view_extension, max_view_extension); } diff --git a/src/buildblock/interpolate_projdata.cxx b/src/buildblock/interpolate_projdata.cxx index 817bfcd3fd..2e9b63ce58 100644 --- a/src/buildblock/interpolate_projdata.cxx +++ b/src/buildblock/interpolate_projdata.cxx @@ -20,9 +20,10 @@ */ #include "stir/ProjData.h" //#include "stir/display.h" +#include "stir/ProjDataInMemory.h" +#include "stir/ProjDataInterfile.h" #include "stir/ProjDataInfo.h" #include "stir/ProjDataInfoCylindricalNoArcCorr.h" -#include "stir/IndexRange.h" #include "stir/BasicCoordinate.h" #include "stir/Sinogram.h" #include "stir/SegmentBySinogram.h" @@ -33,8 +34,11 @@ #include "stir/interpolate_axial_position.h" #include "stir/extend_projdata.h" #include "stir/numerics/sampling_functions.h" +#include "stir/Coordinate4D.h" +#include "stir/IndexRange4D.h" +#include "stir/LORCoordinates.h" #include - +#include START_NAMESPACE_STIR namespace detail_interpolate_projdata @@ -170,7 +174,8 @@ using namespace detail_interpolate_projdata; Succeeded interpolate_projdata(ProjData& proj_data_out, - const ProjData& proj_data_in, const BSpline::BSplineType these_types, + const shared_ptr proj_data_in, + const BSpline::BSplineType &these_types, const bool remove_interleaving, const bool use_view_offset) { @@ -182,7 +187,7 @@ interpolate_projdata(ProjData& proj_data_out, Succeeded interpolate_projdata(ProjData& proj_data_out, - const ProjData& proj_data_in, + const shared_ptr proj_data_in, const BasicCoordinate<3, BSpline::BSplineType> & these_types, const bool remove_interleaving, const bool use_view_offset) @@ -192,7 +197,7 @@ interpolate_projdata(ProjData& proj_data_out, warning("interpolate_projdata with use_view_offset is EXPERIMENTAL and NOT TESTED."); const ProjDataInfo & proj_data_in_info = - *proj_data_in.get_proj_data_info_sptr(); + *proj_data_in->get_proj_data_info_sptr(); const ProjDataInfo & proj_data_out_info = *proj_data_out.get_proj_data_info_sptr(); @@ -204,7 +209,7 @@ interpolate_projdata(ProjData& proj_data_out, if (proj_data_in_info.get_scanner_sptr()->get_scanner_geometry()=="BlocksOnCylindrical") { - interpolate_axial_position(proj_data_out,proj_data_in); + interpolate_axial_position(proj_data_out,*proj_data_in); return Succeeded::yes; } // check for the same ring radius @@ -268,10 +273,10 @@ interpolate_projdata(ProjData& proj_data_out, const SegmentBySinogram non_interleaved_segment = make_non_interleaved_segment(*non_interleaved_proj_data_info_sptr, - proj_data_in.get_segment_by_sinogram(0)); + proj_data_in->get_segment_by_sinogram(0)); // display(non_interleaved_segment, non_interleaved_segment.find_max(),"non-inter"); - Array<3,float> extended = - extend_segment_in_views(non_interleaved_segment, 2, 2); + Array<3,float> extended = + extend_segment_in_views(non_interleaved_segment, 2, 2, non_interleaved_segment); for (int z=extended.get_min_index(); z<= extended.get_max_index(); ++z) { for (int y=extended[z].get_min_index(); y<= extended[z].get_max_index(); ++y) @@ -286,9 +291,9 @@ interpolate_projdata(ProjData& proj_data_out, proj_data_interpolator.set_coef(extended); } else - { - Array<3,float> extended = - extend_segment_in_views(proj_data_in.get_segment_by_sinogram(0), 2, 2); + { //TODO: be removed ... + Array<3,float> extended = + extend_segment_in_views(proj_data_in->get_segment_by_sinogram(0), 2, 2, proj_data_in->get_segment_by_sinogram(0)); for (int z=extended.get_min_index(); z<= extended.get_max_index(); ++z) { for (int y=extended[z].get_min_index(); y<= extended[z].get_max_index(); ++y) @@ -313,4 +318,348 @@ interpolate_projdata(ProjData& proj_data_out, return Succeeded::yes; } +Succeeded +interpolate_projdata_3d(shared_ptr projdata_out, + const shared_ptr projdata_in_sptr, + const BSpline::BSplineType & spline_type, + const bool use_view_offset) +{ + if (use_view_offset) + warning("interpolate_projdata with use_view_offset is EXPERIMENTAL and NOT TESTED."); + + const shared_ptr projdata_in_info_sptr = + projdata_in_sptr->get_proj_data_info_sptr(); + + shared_ptr projdata_in_up_info_sptr(projdata_in_sptr->get_proj_data_info_sptr()->clone()); + projdata_in_up_info_sptr->set_num_views(projdata_out->get_num_views()); + projdata_in_up_info_sptr->set_num_tangential_poss(projdata_out->get_num_tangential_poss()); + + if (projdata_in_info_sptr->get_scanner_sptr()->get_scanner_geometry()=="BlocksOnCylindrical") + { + interpolate_axial_position(*projdata_out, *projdata_in_sptr); + return Succeeded::yes; + } + + // check for the same ring radius + // This is strictly speaking only necessary for non-arccorrected data, but + // we leave it in for all cases. + if (fabs(projdata_in_info_sptr->get_scanner_sptr()->get_inner_ring_radius() - + projdata_out->get_proj_data_info_sptr()->get_scanner_sptr()->get_inner_ring_radius()) > 1) + { + error("interpolate_projdata needs both projection to be of a scanner with the same ring radius"); + } + + std::string output_filename = "tmp_in_up"; + // ProjDataInMemory in_up_projdata(projdata_in_sptr->get_exam_info_sptr(), + // projdata_in_info_sptr->create_shared_clone(), + // 1); // I pressume 1 but we should check! + // For larger scanners ProjDataInterfile might be the right option. + ProjDataInterfile in_up_projdata(projdata_in_sptr->get_exam_info_sptr(), + projdata_in_up_info_sptr, + output_filename, std::ios::out); + + info("interpolate_projdata: Interpolating views and tangential positions ..."); + + bool flag_something_went_wrong = false; + + #ifdef STIR_OPENMP + # if _OPENMP <201107 + #pragma omp parallel for + # else + #pragma omp parallel for schedule(dynamic) + # endif + #endif +#ifdef STIR_TOF + for(int i_tof_in = projdata_in_sptr->get_min_tof_pos_num(); + i_tof_in <= projdata_in_sptr->get_max_tof_pos_num(); ++i_tof_in) + { +#endif + for(int i_seg_in = projdata_in_sptr->get_min_segment_num(); + i_seg_in <= projdata_in_sptr->get_max_segment_num(); ++i_seg_in) + { + info(boost::format("Now processing segment #: %1%") % i_seg_in); + + BasicCoordinate<2, BSpline::BSplineType> these_types2; + these_types2[1]=these_types2[2]=spline_type; + BasicCoordinate<2, double> offset, step; + + const float in_sampling_phi = + (projdata_in_info_sptr->get_phi(Bin(i_seg_in,1,0,0/*, i_tof_in*/)) - projdata_in_info_sptr->get_phi(Bin(i_seg_in,0,0,0/*, i_tof_in*/))); + const float in_up_sampling_phi = + projdata_out->get_proj_data_info_sptr()->get_phi(Bin(i_seg_in,1,0,0/*, i_tof_in*/)) - projdata_out->get_proj_data_info_sptr()->get_phi(Bin(i_seg_in,0,0,0/*, i_tof_in*/)); + + const float in_view_offset = 0; +// use_view_offset +// ? projdata_in_info_sptr->get_scanner_ptr()->get_intrinsic_azimuthal_tilt() +// : 0.F; + const float in_up_view_offset = 0; +// use_view_offset +// ? projdata_in_up_info_sptr->get_scanner_ptr()->get_intrinsic_azimuthal_tilt() +// : 0.F; + + offset[1] = + (projdata_in_info_sptr->get_phi(Bin(i_seg_in,0,0,0/*, i_tof_in*/)) + in_view_offset - + projdata_out->get_proj_data_info_sptr()->get_phi(Bin(i_seg_in,0,0,0/*, i_tof_in*/)) - in_up_view_offset) / + in_sampling_phi; + step[1] = + in_up_sampling_phi/in_sampling_phi; + + const float in_sampling_s = projdata_in_info_sptr->get_sampling_in_s(Bin(i_seg_in,0,0,0/*, i_tof_in*/)); + const float in_up_sampling_s = projdata_out->get_proj_data_info_sptr()->get_sampling_in_s(Bin(i_seg_in,0,0,0/*, i_tof_in*/)); + + float dd = projdata_out->get_proj_data_info_sptr()->get_s(Bin(i_seg_in,0, 0,0/*, i_tof_in*/)); + float ee = projdata_in_info_sptr->get_s(Bin(i_seg_in,0,0, 0/*, i_tof_in*/)); + + offset[2] = + ( dd-ee ) / in_sampling_s; + step[2]= + in_up_sampling_s / in_sampling_s; + + std::cout << "n3" << std::endl; + int inv_seg = -i_seg_in; + + Array<3,float> extended3 = + extend_segment_in_views(projdata_in_sptr->get_segment_by_sinogram(i_seg_in/*, i_tof_in*/), 2, 2, + projdata_in_sptr->get_segment_by_sinogram(inv_seg/*, i_tof_in*/)); + SegmentBySinogram sino_3D_out = in_up_projdata.get_proj_data_info_sptr()->get_empty_segment_by_sinogram(i_seg_in); + + std::cout << "n6" << std::endl; + for(int i_axial = projdata_in_sptr->get_min_axial_pos_num(i_seg_in); + i_axial <= projdata_in_sptr->get_max_axial_pos_num(i_seg_in); ++i_axial) + { + Array<2,float> extended = extended3[i_axial]; + BSpline::BSplinesRegularGrid<2, float, float> proj_data_interpolator(these_types2); + for (int y=extended.get_min_index(); y<= extended.get_max_index(); ++y) + { + const int old_min = extended[y].get_min_index(); + const int old_max = extended[y].get_max_index(); + extended[y].grow(old_min-1, old_max+1); + extended[y][old_min-1] = extended[y][old_min]; + extended[y][old_max+1] = extended[y][old_max]; + } + std::cout << "n7" << std::endl; + proj_data_interpolator.set_coef(extended); + std::cout << "n7.2" << std::endl; + Sinogram sino_2D_out = sino_3D_out.get_sinogram(i_axial); + sample_function_on_regular_grid(sino_2D_out, proj_data_interpolator, offset, step); + sino_3D_out.set_sinogram(sino_2D_out, i_axial); + + } + if (in_up_projdata.set_segment(sino_3D_out) == Succeeded::no) + { + flag_something_went_wrong = true; + } + + } + +#ifdef STIR_TOF + } +#endif + + if(flag_something_went_wrong == true) + error("interpolate_projdata: Something went wrong in the first level of interpolation"); + + info("interpolate_projdata: Finished interpolating views and tangential positions!"); + + const stir::shared_ptr< ProjData > in_up_projdata_in = + ProjData::read_from_file("tmp_in_up.hs"); + // std::cout << in_up_projdata_in->get_proj_data_info_sptr()->parameter_info() << std::endl; + + + // std::cout << projdata_out->get_proj_data_info_sptr()->parameter_info() << std::endl; + +// if(projdata_out->get_num_segments() == 1) +// { +// SegmentByView v = in_up_projdata_in->get_segment_by_view(0); +// projdata_out->set_segment(v); +// int summm = v.sum(); + +// std::cout << "in_up_projdata_in sum: " << summm << std::endl; +// summm= projdata_out->get_segment_by_sinogram(0).sum(); +// std::cout << "projdata_out sum: " << summm << std::endl; +// return Succeeded::yes; +// } + + + const shared_ptr in_up_projdata_in_info_sptr = + dynamic_pointer_cast(in_up_projdata_in->get_proj_data_info_sptr()->create_shared_clone()); + + if(is_null_ptr(in_up_projdata_in_info_sptr)) + error("BlockGeometry or arc correction, whye are we here??"); + + // Compressed data are not supported + if (in_up_projdata_in_info_sptr->get_min_ring_difference(0) != in_up_projdata_in_info_sptr->get_max_ring_difference(0)) + { + warning("Compressed data are not supported, yet!"); + return Succeeded::no; + } + + + const shared_ptr out_projdata_info_sptr = + dynamic_pointer_cast(projdata_out->get_proj_data_info_sptr()->create_shared_clone()); + + if(is_null_ptr(out_projdata_info_sptr)) + error("Interpolate_projdata: Interpolation in 3D needs non arc corrected data."); + + // const stir::shared_ptr projdata_out_info = projdata_out->get_proj_data_info_sptr(); + + // Add two more rings, to have a wider space for the interpolation + IndexRange4D mich_index( 0, in_up_projdata_in_info_sptr->get_scanner_sptr()->get_num_rings()+1, + 0, in_up_projdata_in_info_sptr->get_scanner_sptr()->get_num_rings()+1, + 0, out_projdata_info_sptr->get_scanner_sptr()->get_num_detectors_per_ring()-1, + 0, out_projdata_info_sptr->get_scanner_sptr()->get_num_detectors_per_ring()-1); + + BasicCoordinate<4, BSpline::BSplineType> these_types_4; + these_types_4[1]=these_types_4[2]=these_types_4[3]=these_types_4[4] = spline_type; + +#ifdef STIR_TOF // Now it is better to NOT parallelise the TOF bins but the geometric loops. + info("interpolate_projdata: Creating michelogram for 3D interpolation ..."); + for(int i_tof_in = in_up_projdata_in->get_min_tof_pos_num(); + i_tof_in <= in_up_projdata_in->get_max_tof_pos_num(); ++i_tof_in) + { + const int cur_tof = i_tof_in; +#endif + + // Create the 4D Michelogram R1.R1.D1.D2 + Array<4, float> downsampled_array_4d(mich_index); + downsampled_array_4d.fill(0.0); + +#ifdef STIR_OPENMP +#pragma omp parallel for schedule(dynamic) //collapse(3) +#endif + for (int i_seg = in_up_projdata_in_info_sptr->get_min_segment_num(); + i_seg <= in_up_projdata_in_info_sptr->get_max_segment_num(); ++i_seg) + { + int min_axial_pos = in_up_projdata_in_info_sptr->get_min_axial_pos_num(i_seg); + int max_axial_pos = in_up_projdata_in_info_sptr->get_max_axial_pos_num(i_seg); + + // info( std::string("interpolate_projdata: Processing segment:") << i_seg ); + std::cout <<"interpolate_projdata: Processing segment:" << i_seg << std::endl; + + SegmentBySinogram sino3D = in_up_projdata_in->get_segment_by_sinogram(i_seg/*, cur_tof*/); + + for(int i_axial = min_axial_pos; i_axial <= max_axial_pos; ++i_axial) + { + int r1, r2; + in_up_projdata_in_info_sptr->get_ring_pair_for_segment_axial_pos_num(r1, r2, i_seg, i_axial); + // For the view and tangential position we have to use the finer template, as otherwise we might run + // into problems with the downsampled. + for(int i_view = in_up_projdata_in_info_sptr->get_min_view_num(); i_view <= in_up_projdata_in_info_sptr->get_max_view_num(); ++i_view) + { + for (int i_tang = in_up_projdata_in_info_sptr->get_min_tangential_pos_num(); i_tang <= in_up_projdata_in_info_sptr->get_max_tangential_pos_num(); ++i_tang) + { + int c1,c2; + out_projdata_info_sptr->get_det_num_pair_for_view_tangential_pos_num(c1, c2, i_view, i_tang); + // dynamic_cast(in_up_projdata.get_proj_data_info_sptr().get())->get_det_pair_for_bin(c1, r1, c2, r2, tmp_bin); + // std::cout << r2 << " " << r1 << " " << c2 << " " << c1 << std::endl; + downsampled_array_4d[r1+1][r2+1][c1][c2] = sino3D[i_axial][i_view][i_tang]; + } + } + } + } + std::cout << "kkd" << std::endl; + // Fill the extra rings with the values from their neigbours + int last_ring = in_up_projdata_in_info_sptr->get_scanner_sptr()->get_num_rings(); + for(int i_c1 = 0; i_c1 < out_projdata_info_sptr->get_scanner_sptr()->get_num_detectors_per_ring(); ++i_c1) + { + for(int i_c2 = 0; i_c2 < out_projdata_info_sptr->get_scanner_sptr()->get_num_detectors_per_ring(); ++i_c2) + { + downsampled_array_4d[0][0][i_c1][i_c2] = + downsampled_array_4d[1][1][i_c1][i_c2]; + downsampled_array_4d[last_ring+1][last_ring+1][i_c1][i_c2] = + downsampled_array_4d[last_ring][last_ring][i_c1][i_c2]; + } + + } + + BSpline::BSplinesRegularGrid<4, float, float> mich_data_interpolator(spline_type); + mich_data_interpolator.set_coef(downsampled_array_4d); + + info("interpolate_projdata: Finished michelogram!"); + + const float angle_spacing = 1/(in_up_projdata_in_info_sptr->get_phi(Bin(0,1,0, 0/*, cur_tof*/))- + in_up_projdata_in_info_sptr->get_phi(Bin(0,0,0,0/*, cur_tof*/))); + const float eff_ring_radius = out_projdata_info_sptr->get_scanner_sptr()->get_effective_ring_radius(); + +#ifdef STIR_OPENMP +#pragma omp parallel for schedule(dynamic) +#endif + for (int i_seg = out_projdata_info_sptr->get_min_segment_num(); i_seg <= out_projdata_info_sptr->get_max_segment_num(); ++i_seg) + { + int min_axial_pos = out_projdata_info_sptr->get_min_axial_pos_num(i_seg); + int max_axial_pos = out_projdata_info_sptr->get_max_axial_pos_num(i_seg); + SegmentBySinogram _sino3D = out_projdata_info_sptr->get_empty_segment_by_sinogram(i_seg, false/*, cur_tof*/); + //info(boost::format("interpolate_projdata: Interpolating michelogram, segment: %1%") % i_seg); + + for(int i_axial = min_axial_pos; i_axial <= max_axial_pos; ++i_axial) + { + const double ring_spacing_in = 1/in_up_projdata_in_info_sptr->get_ring_spacing(); + const double ring_spacing_out = out_projdata_info_sptr->get_ring_spacing(); + const float offset = 0.5*in_up_projdata_in_info_sptr->get_ring_spacing() + out_projdata_info_sptr->get_scanner_sptr()->get_num_rings() * ring_spacing_out * 0.5; + + for(int i_view = out_projdata_info_sptr->get_min_view_num(); i_view <= out_projdata_info_sptr->get_max_view_num(); ++i_view) + { + for (int i_tang = out_projdata_info_sptr->get_min_tangential_pos_num(); i_tang <= out_projdata_info_sptr->get_max_tangential_pos_num(); ++i_tang) + { + //std::cout << i_axial << " " << i_view << " " << i_tang << std::endl; + float value = 0.0; + Bin tmp_bin(i_seg, i_view, i_axial, i_tang, /*cur_tof,*/ value); + + LORInAxialAndNoArcCorrSinogramCoordinates lor_sino; + LORInCylinderCoordinates lor_cyl; + + out_projdata_info_sptr->get_LOR(lor_sino, tmp_bin); + int d1, d2, dr1, dr2; + out_projdata_info_sptr->get_det_pair_for_bin(d1, dr1, d2, dr2, tmp_bin); + lor_sino.change_representation(lor_cyl, eff_ring_radius); + + double c1(lor_cyl.p1().psi()*angle_spacing); + double c2(lor_cyl.p2().psi()*angle_spacing); + double r1((lor_cyl.p1().z() + offset )*(ring_spacing_in)); + double r2((lor_cyl.p2().z() + offset )*(ring_spacing_in)); + + //int cur_tilt, cur_slice, cur_phi, cur_rad; + //sinoMap.getSinoCoordsForCrystals( + // 0, d1, dr1, + // 0, d2, dr2, + // &cur_tilt, &cur_slice, &cur_phi, &cur_rad, 0); + + + // if(!is_gap(dr1) && !is_gap(dr2)) + { + // std::cout << r1 << " " << r2 << " " << c1 << " " << c2 << std::endl; + BasicCoordinate<4, double> pos = make_coordinate(r1, r2, c1, c2); + //if(cur_phi * mh_normsino.numray + cur_rad > sliceSize || + //cur_slice > norm_seg.size()) + //std::cout << "Nikos: Overflow!" << cur_slice << " of " << norm_seg.size() << ", " << cur_phi*mh_normsino.numray + cur_rad <<" of " << sliceSize << std::endl; + value = mich_data_interpolator(pos);// * norm_seg[cur_slice][cur_phi * mh_normsino.numray + cur_rad]; + } + + _sino3D[i_axial][i_view][i_tang] = value; + } + } + } + projdata_out->set_segment(_sino3D); + } + +#ifdef STIR_TOF + } +#endif + + info("interpolate_projdata: Finished michelogram interpolation!"); + return Succeeded::yes; +} + + + + + + + + + + + + + END_NAMESPACE_STIR diff --git a/src/buildblock/scale_sinograms.cxx b/src/buildblock/scale_sinograms.cxx index fc52e460e0..c51b02a8ce 100644 --- a/src/buildblock/scale_sinograms.cxx +++ b/src/buildblock/scale_sinograms.cxx @@ -59,7 +59,7 @@ scale_sinograms( Array<2,float> -get_scale_factors_per_sinogram(const ProjData& numerator_proj_data, +get_scale_factors_per_sinogram(const ProjData& numerator_proj_data, const ProjData& denominator_proj_data, const ProjData& weights_proj_data) { @@ -116,7 +116,7 @@ get_scale_factors_per_sinogram(const ProjData& numerator_proj_data, { warning("Problem at segment %d, axial pos %d in finding sinogram scaling factor.\n" "Weighted data in denominator %g is very small compared to total in sinogram %g.\n" - "Adjust weights?.\n" + "Adjust weights?. This could be an axial gap.\n" "I will use scale factor %g", bin.segment_num(),bin.axial_pos_num(), total_in_denominator[bin.segment_num()][bin.axial_pos_num()], diff --git a/src/include/stir/extend_projdata.h b/src/include/stir/extend_projdata.h index 53ef731886..f07d510582 100644 --- a/src/include/stir/extend_projdata.h +++ b/src/include/stir/extend_projdata.h @@ -31,10 +31,12 @@ This is probably only useful before calling interpolation routines, or for FORE. */ Array<3,float> extend_segment_in_views(const SegmentBySinogram& sino, - const int min_view_extension, const int max_view_extension); + const int min_view_extension, const int max_view_extension, + const SegmentBySinogram& inv_sino); Array<2,float> extend_sinogram_in_views(const Sinogram& sino, - const int min_view_extension, const int max_view_extension); + const int min_view_extension, const int max_view_extension, + const Sinogram& inv_sino); //@} END_NAMESPACE_STIR diff --git a/src/include/stir/interpolate_projdata.h b/src/include/stir/interpolate_projdata.h index 114f533a2c..a8d1831bea 100644 --- a/src/include/stir/interpolate_projdata.h +++ b/src/include/stir/interpolate_projdata.h @@ -22,6 +22,7 @@ START_NAMESPACE_STIR class ProjData; +class ProjDataInfo; template class BasicCoordinate; template class Sinogram; template class SegmentBySinogram; @@ -51,16 +52,22 @@ template class SegmentBySinogram; //@{ Succeeded interpolate_projdata(ProjData& proj_data_out, - const ProjData& proj_data_in, - const BSpline::BSplineType spline_type, + const shared_ptr proj_data_in, + const BSpline::BSplineType& spline_type, const bool remove_interleaving = false, const bool use_view_offset = false); Succeeded interpolate_projdata(ProjData& proj_data_out, - const ProjData& proj_data_in, + const shared_ptr proj_data_in, const BasicCoordinate<3, BSpline::BSplineType> & spline_type, const bool remove_interleaving = false, - const bool use_view_offset = false); + const bool use_view_offset = false); + +Succeeded +interpolate_projdata_3d(shared_ptr projdata_out, + const shared_ptr projdata_in_sptr, + const BSpline::BSplineType & spline_type, + const bool use_view_offset = false); //@} END_NAMESPACE_STIR diff --git a/src/include/stir/numerics/sampling_functions.h b/src/include/stir/numerics/sampling_functions.h index ddc2f75668..b112f75cf9 100644 --- a/src/include/stir/numerics/sampling_functions.h +++ b/src/include/stir/numerics/sampling_functions.h @@ -18,6 +18,8 @@ */ +#include "stir/Array.h" +#include "stir/common.h" START_NAMESPACE_STIR /*! @@ -49,6 +51,13 @@ void sample_function_on_regular_grid(Array<3,elemT>& out, const BasicCoordinate<3, positionT>& offset, const BasicCoordinate<3, positionT>& step); +template +inline +void sample_function_on_regular_grid(Array<2,elemT>& out, + FunctionType func, + const BasicCoordinate<2, positionT>& offset, + const BasicCoordinate<2, positionT>& step); + END_NAMESPACE_STIR #include "stir/numerics/sampling_functions.inl" diff --git a/src/include/stir/numerics/sampling_functions.inl b/src/include/stir/numerics/sampling_functions.inl index b504727f97..d1635b64fe 100644 --- a/src/include/stir/numerics/sampling_functions.inl +++ b/src/include/stir/numerics/sampling_functions.inl @@ -32,26 +32,64 @@ void sample_function_on_regular_grid(Array<3,elemT>& out, BasicCoordinate<3, positionT> relative_positions; index_out[1]=min_out[1]; relative_positions[1]= index_out[1] * step[1] - offset[1] ; - const BasicCoordinate<3, positionT> max_relative_positions= - (BasicCoordinate<3,positionT>(max_out)+static_cast(.001)) * step + offset; + const BasicCoordinate<3, positionT> max_relative_positions= + (BasicCoordinate<3,positionT>(max_out)+static_cast(.001)) * step + offset; for (; index_out[1]<=max_out[1] && relative_positions[1]<=max_relative_positions[1]; ++index_out[1], relative_positions[1]+= step[1]) - { + { index_out[2]=min_out[2]; - relative_positions[2]= index_out[2] * step[2] + offset[2] ; + relative_positions[2]= index_out[2] * step[2] + offset[2] ; for (; index_out[2]<=max_out[2] && relative_positions[2]<=max_relative_positions[2]; ++index_out[2], relative_positions[2]+= step[2]) - { + { index_out[3]=min_out[3]; - relative_positions[3]= index_out[3] * step[3] + offset[3] ; + relative_positions[3]= index_out[3] * step[3] + offset[3] ; for (; index_out[3]<=max_out[3] && relative_positions[3]<=max_relative_positions[3]; - ++index_out[3], relative_positions[3]+= step[3]) - out[index_out] = func(relative_positions) ; - } - } + ++index_out[3], relative_positions[3]+= step[3]) + { + out[index_out] = func(relative_positions) ; + } + } + } +} + +template +void sample_function_on_regular_grid(Array<2,elemT>& out, + FunctionType func, + const BasicCoordinate<2, positionT>& offset, + const BasicCoordinate<2, positionT>& step) +{ + BasicCoordinate<2,int> min_out, max_out; + IndexRange<2> out_range = out.get_index_range(); + + if (!out_range.get_regular_range(min_out,max_out)) + warning("Output must be regular range!"); + +// in_total = out.sum(); + BasicCoordinate<2, int> index_out; + BasicCoordinate<2, positionT> relative_positions; + index_out[1] = min_out[1]; + relative_positions[1] = index_out[1] * step[1] - offset[1] ; + + const BasicCoordinate<2, positionT> max_relative_positions= + (BasicCoordinate<2,positionT>(max_out)+static_cast(.001)) * step + offset; + for (; + index_out[1]<=max_out[1] && relative_positions[1]<=max_relative_positions[1]; + ++index_out[1], relative_positions[1]+= step[1]) + { + index_out[2]=min_out[2]; + relative_positions[2]= index_out[2] * step[2] + offset[2] ; + for (; + index_out[2]<=max_out[2] && relative_positions[2]<=max_relative_positions[2]; + ++index_out[2], relative_positions[2]+= step[2]) + { + + out[index_out] = func(relative_positions) ; + } + } } END_NAMESPACE_STIR diff --git a/src/include/stir/scatter/ScatterEstimation.h b/src/include/stir/scatter/ScatterEstimation.h index 3e354ec435..da5dc22cf1 100644 --- a/src/include/stir/scatter/ScatterEstimation.h +++ b/src/include/stir/scatter/ScatterEstimation.h @@ -97,14 +97,15 @@ class ScatterEstimation: public ParsingObject static void upsample_and_fit_scatter_estimate(ProjData& scaled_scatter_proj_data, const ProjData& emission_proj_data, - const ProjData& scatter_proj_data, + shared_ptr scatter_proj_data, BinNormalisation& scatter_normalisation, const ProjData& weights_proj_data, const float min_scale_factor, const float max_scale_factor, const unsigned half_filter_width, BSpline::BSplineType spline_type = BSpline::BSplineType::linear, - const bool remove_interleaving = true); + const bool remove_interleaving = true, + const bool do_3D = false); //! Default constructor (calls set_defaults()) diff --git a/src/scatter_buildblock/ScatterEstimation.cxx b/src/scatter_buildblock/ScatterEstimation.cxx index a86b96b94c..a72352f6a1 100644 --- a/src/scatter_buildblock/ScatterEstimation.cxx +++ b/src/scatter_buildblock/ScatterEstimation.cxx @@ -473,7 +473,9 @@ set_up() error("ScatterEstimation: Please define a scatter simulation method. Aborting."); if (!run_in_2d_projdata) - error("ScatterEstimation: Currently, only running the estimation in 2D is supported."); + warning("ScatterEstimation: Currently, only running the estimation in 2D is supported."); + else + warning("This is highly experimental. Try it on your own risk!"); if(!this->recompute_mask_projdata) { @@ -496,16 +498,18 @@ set_up() #if 1 // Calculate the SSRB - if (input_projdata_sptr->get_num_segments() > 1) - { - info("ScatterEstimation: Running SSRB on input data..."); - this->input_projdata_2d_sptr = make_2D_projdata_sptr(this->input_projdata_sptr); - } - else - { - input_projdata_2d_sptr = input_projdata_sptr; - } - + if (run_in_2d_projdata) + { + if (input_projdata_sptr->get_num_segments() > 1) + { + info("ScatterEstimation: Running SSRB on input data..."); + this->input_projdata_2d_sptr = make_2D_projdata_sptr(this->input_projdata_sptr); + } + else + { + input_projdata_2d_sptr = input_projdata_sptr; + } + } #else { std::string tmp_input2D = "./extras/nema_proj_f1g1d0b0.hs_2d.hs"; @@ -642,92 +646,113 @@ set_up_iterative(shared_ptr iterative_object->set_input_data(this->input_projdata_sptr); - - // - // Multiplicative projdata - // - shared_ptr tmp_atten_projdata_sptr = - this->get_attenuation_correction_factors_sptr(this->multiplicative_binnorm_sptr); - shared_ptr atten_projdata_2d_sptr; - - info("ScatterEstimation: 3.Calculating the attenuation projection data..."); - - if( tmp_atten_projdata_sptr->get_num_segments() > 1) + if (run_in_2d_projdata) { - info("ScatterEstimation: Running SSRB on attenuation correction coefficients ..."); + info("ScatterEstimation: 3.Calculating the attenuation projection data..."); + // + // Multiplicative projdata + // + shared_ptr tmp_atten_projdata_sptr = + this->get_attenuation_correction_factors_sptr(this->multiplicative_binnorm_sptr); + shared_ptr atten_projdata_2d_sptr; + + if( tmp_atten_projdata_sptr->get_num_segments() > 1) + { + info("ScatterEstimation: 3.1. Running SSRB on attenuation correction coefficients ..."); - std::string out_filename = "tmp_atten_sino_2d.hs"; - atten_projdata_2d_sptr=make_2D_projdata_sptr(tmp_atten_projdata_sptr, out_filename); - - } - else - { - // TODO: this needs more work. -- Setting directly 2D proj_data is buggy right now. - atten_projdata_2d_sptr = tmp_atten_projdata_sptr; - } + std::string out_filename = "tmp_atten_sino_2d.hs"; + atten_projdata_2d_sptr=make_2D_projdata_sptr(tmp_atten_projdata_sptr, out_filename); - info("ScatterEstimation: 4.Calculating the normalisation data..."); - { - if (run_in_2d_projdata) + } + else { - shared_ptr norm3d_sptr = - this->get_normalisation_object_sptr(this->multiplicative_binnorm_sptr); - shared_ptr norm_coeff_2d_sptr; + // TODO: this needs more work. -- Setting directly 2D proj_data is buggy right now. + atten_projdata_2d_sptr = tmp_atten_projdata_sptr; + } + info("ScatterEstimation: 4.Calculating the normalisation data..."); + shared_ptr norm3d_sptr = + this->get_normalisation_object_sptr(this->multiplicative_binnorm_sptr); + shared_ptr norm_coeff_2d_sptr; - if ( input_projdata_sptr->get_num_segments() > 1) - { - // Some BinNormalisation classes don't know about SSRB. - // we need to get norm2d=1/SSRB(1/norm3d)) + if ( input_projdata_sptr->get_num_segments() > 1) + { + // Some BinNormalisation classes don't know about SSRB. + // we need to get norm2d=1/SSRB(1/norm3d)) - info("ScatterEstimation: Constructing 2D normalisation coefficients ..."); + info("ScatterEstimation: 4.1 Constructing 2D normalisation coefficients ..."); - std::string out_filename = "tmp_inverted_normdata.hs"; - shared_ptr inv_projdata_3d_sptr = create_new_proj_data(out_filename, - this->input_projdata_sptr->get_exam_info_sptr(), - this->input_projdata_sptr->get_proj_data_info_sptr()->create_shared_clone()); - inv_projdata_3d_sptr->fill(1.f); + std::string out_filename = "tmp_inverted_normdata.hs"; + shared_ptr inv_projdata_3d_sptr = create_new_proj_data(out_filename, + this->input_projdata_sptr->get_exam_info_sptr(), + this->input_projdata_sptr->get_proj_data_info_sptr()->create_shared_clone()); + inv_projdata_3d_sptr->fill(1.f); - out_filename = "tmp_normdata_2d.hs"; - shared_ptr norm_projdata_2d_sptr = create_new_proj_data(out_filename, - this->input_projdata_2d_sptr->get_exam_info_sptr(), - this->input_projdata_2d_sptr->get_proj_data_info_sptr()->create_shared_clone()); - norm_projdata_2d_sptr->fill(0.f); + // Essentially since inv_projData_sptr is 1s then this is an inversion. + // inv_projdata_sptr = 1/norm3d + norm3d_sptr->undo(*inv_projdata_3d_sptr); - // Essentially since inv_projData_sptr is 1s then this is an inversion. - // inv_projdata_sptr = 1/norm3d - norm3d_sptr->undo(*inv_projdata_3d_sptr); + info("ScatterEstimation: 4.2 Performing SSRB on efficiency factors ..."); - info("ScatterEstimation: Performing SSRB on efficiency factors ..."); + shared_ptr norm_projdata_2d_sptr=make_2D_projdata_sptr(inv_projdata_3d_sptr); + norm_projdata_2d_sptr->write_to_file("extras/tmp_normdata_2d.hs"); - norm_projdata_2d_sptr=make_2D_projdata_sptr(inv_projdata_3d_sptr); + // Crucial: Avoid divisions by zero!! + // This should be resolved after https://github.com/UCL/STIR/issues/348 + pow_times_add min_threshold (0.0f, 1.0f, 1.0f, 1E-20f, NumericInfo().max_value()); + apply_to_proj_data(*norm_projdata_2d_sptr, min_threshold); - // Crucial: Avoid divisions by zero!! - // This should be resolved after https://github.com/UCL/STIR/issues/348 - pow_times_add min_threshold (0.0f, 1.0f, 1.0f, 1E-20f, NumericInfo().max_value()); - apply_to_proj_data(*norm_projdata_2d_sptr, min_threshold); + pow_times_add invert (0.0f, 1.0f, -1.0f, NumericInfo().min_value(), NumericInfo().max_value()); + apply_to_proj_data(*norm_projdata_2d_sptr, invert); - pow_times_add invert (0.0f, 1.0f, -1.0f, NumericInfo().min_value(), NumericInfo().max_value()); - apply_to_proj_data(*norm_projdata_2d_sptr, invert); + norm_coeff_2d_sptr.reset(new BinNormalisationFromProjData(norm_projdata_2d_sptr)); + } + else + { + norm_coeff_2d_sptr = norm3d_sptr; + } - norm_coeff_2d_sptr.reset(new BinNormalisationFromProjData(norm_projdata_2d_sptr)); - } - else - { - norm_coeff_2d_sptr = norm3d_sptr; - } + shared_ptratten_coeff_2d_sptr(new BinNormalisationFromProjData(atten_projdata_2d_sptr)); + this->multiplicative_binnorm_2d_sptr.reset( + new ChainedBinNormalisation(norm_coeff_2d_sptr, atten_coeff_2d_sptr)); - shared_ptratten_coeff_2d_sptr(new BinNormalisationFromProjData(atten_projdata_2d_sptr)); - this->multiplicative_binnorm_2d_sptr.reset( - new ChainedBinNormalisation(norm_coeff_2d_sptr, atten_coeff_2d_sptr)); + this->multiplicative_binnorm_2d_sptr->set_up(this->back_projdata_sptr->get_exam_info_sptr(), this->input_projdata_2d_sptr->get_proj_data_info_sptr()->create_shared_clone()); + iterative_object->get_objective_function_sptr()->set_normalisation_sptr(multiplicative_binnorm_2d_sptr); + } + else + { +// shared_ptr norm_coeff_3d_sptr; +// // run_in_2d_projdata +// // I need the inverted norm factors here!!!! << Tomorrow +// // Create the inverted version of the normalisation factors with very large values in the gaps. +// //As for reconstruction. +// shared_ptr norm3d_sptr = +// this->get_normalisation_object_sptr(this->multiplicative_binnorm_sptr); +// std::string out_filename = "tmp_inverted_normdata.hs"; +// shared_ptr inv_projdata_3d_sptr = create_new_proj_data(out_filename, +// this->input_projdata_sptr->get_exam_info_sptr(), +// this->input_projdata_sptr->get_proj_data_info_sptr()->create_shared_clone()); +// inv_projdata_3d_sptr->fill(1.f); + +// norm3d_sptr->undo(*inv_projdata_3d_sptr); + +// // Crucial: Avoid divisions by zero!! +// // This should be resolved after https://github.com/UCL/STIR/issues/348 +// pow_times_add min_threshold (0.0f, 1.0f, 1.0f, 1E-20f, NumericInfo().max_value()); +// apply_to_proj_data(*inv_projdata_3d_sptr, min_threshold); + +// pow_times_add invert (0.0f, 1.0f, -1.0f, NumericInfo().min_value(), NumericInfo().max_value()); +// apply_to_proj_data(*inv_projdata_3d_sptr, invert); + +// norm_coeff_3d_sptr.reset(new BinNormalisationFromProjData(inv_projdata_3d_sptr)); + +// shared_ptratten_coeff_sptr(new BinNormalisationFromProjData(tmp_atten_projdata_sptr)); +// this->multiplicative_binnorm_3d_sptr.reset( +// new ChainedBinNormalisation(norm_coeff_3d_sptr, atten_coeff_sptr)); +// this->multiplicative_binnorm_3d_sptr->set_up(this->input_projdata_sptr->get_exam_info_sptr(), this->input_projdata_sptr->get_proj_data_info_sptr()->create_shared_clone()); + iterative_object->get_objective_function_sptr()->set_normalisation_sptr(multiplicative_binnorm_sptr); - this->multiplicative_binnorm_2d_sptr->set_up(this->back_projdata_sptr->get_exam_info_sptr(), this->input_projdata_2d_sptr->get_proj_data_info_sptr()->create_shared_clone()); - iterative_object->get_objective_function_sptr()->set_normalisation_sptr(multiplicative_binnorm_2d_sptr); - } - else // run_in_2d_projdata - iterative_object->get_objective_function_sptr()->set_normalisation_sptr(multiplicative_binnorm_sptr); } info("ScatterEstimation: Done normalisation coefficients."); - // // Set background (randoms) projdata // @@ -735,16 +760,13 @@ set_up_iterative(shared_ptr if (!is_null_ptr(this->back_projdata_sptr)) { - if( back_projdata_sptr->get_num_segments() > 1) - { - info("ScatterEstimation: Running SSRB on the background data ..."); - - - this->back_projdata_2d_sptr=make_2D_projdata_sptr(back_projdata_sptr); - } - else + if(run_in_2d_projdata) { - this->back_projdata_2d_sptr = back_projdata_sptr; + if( back_projdata_sptr->get_num_segments() > 1) + { + info("ScatterEstimation: 5.2 Running SSRB on the background data ..."); + this->back_projdata_2d_sptr=make_2D_projdata_sptr(back_projdata_sptr); + } } } else // We will need a background for the scatter, so let's create a simple empty ProjData @@ -769,8 +791,6 @@ set_up_iterative(shared_ptr } } - - if (run_in_2d_projdata) { // Normalise in order to get the additive component @@ -851,7 +871,7 @@ process_data() shared_ptr normalisation_factors_sptr = this->get_normalisation_object_sptr(run_in_2d_projdata ? this->multiplicative_binnorm_2d_sptr - : this->multiplicative_binnorm_sptr); + : this->multiplicative_binnorm_sptr); if(run_in_2d_projdata) { scaled_est_projdata_sptr.reset(new ProjDataInMemory(this->input_projdata_2d_sptr->get_exam_info_sptr(), @@ -881,6 +901,7 @@ process_data() else reconstruct_analytic(0, this->current_activity_image_sptr); +// current_activity_image_sptr->read_from_file("extras/initial_activity_image.hv"); if ( run_debug_mode ) { std::string out_filename = extras_path.get_path() + "initial_activity_image"; @@ -956,12 +977,14 @@ process_data() scaled_est_projdata_sptr->fill(0.F); - upsample_and_fit_scatter_estimate(*scaled_est_projdata_sptr, *data_to_fit_projdata_sptr, - *unscaled_est_projdata_sptr, + upsample_and_fit_scatter_estimate(*scaled_est_projdata_sptr, //output: Now is normalized + *data_to_fit_projdata_sptr, // emission + unscaled_est_projdata_sptr, // unscaled: OK *normalisation_factors_sptr, - *this->mask_projdata_sptr, local_min_scale_value, + *this->mask_projdata_sptr, // weights: OK + local_min_scale_value, local_max_scale_value, this->half_filter_width, - spline_type, true); + spline_type, run_in_2d_projdata, !run_in_2d_projdata); if(this->run_debug_mode) { @@ -973,16 +996,16 @@ process_data() } - // When saving we need to go 3D. + if (this->export_scatter_estimates_of_each_iteration || i_scat_iter == this->num_scatter_iterations ) { shared_ptr temp_scatter_projdata; - + // When saving we need to go 3D. if(run_in_2d_projdata) { - info("ScatterEstimation: upsampling scatter to 3D"); + info("ScatterEstimation: upsampling scatter to 3D"); //this is complicated as the 2d scatter estimate was //"unnormalised" (divided by norm2d), so we need to undo this 2D norm, and put a 3D norm in. //unfortunately, currently the values in the gaps in the @@ -1001,14 +1024,14 @@ process_data() pow_times_add add_scalar (-1e-9f, 1.0f, 1.0f, NumericInfo().min_value(), NumericInfo().max_value()); apply_to_proj_data(*temp_projdata, min_threshold); apply_to_proj_data(*temp_projdata, add_scalar); - // threshold back to 0 to avoid getting tiny negatives (due to numerical precision errors) + // threshold back to 0 to avoid getting tiny negatives (due to numerical precision errors) pow_times_add min_threshold_zero (0.0f, 1.0f, 1.0f, 0.f, NumericInfo().max_value()); apply_to_proj_data(*temp_projdata, min_threshold_zero); // ok, we can multiply with the norm normalisation_factors_sptr->apply(*temp_projdata); - // Create proj_data to save the 3d scatter estimate + // Create proj_data to save the 3d scatter estimate if(!this->output_scatter_estimate_prefix.empty()) { std::stringstream convert; @@ -1023,55 +1046,56 @@ process_data() } else { - // TODO should check if we have one already from previous iteration + // TODO should check if we have one already from previous iteration scatter_estimate_sptr.reset( new ProjDataInMemory(this->input_projdata_sptr->get_exam_info_sptr(), - this->input_projdata_sptr->get_proj_data_info_sptr())); + this->input_projdata_sptr->get_proj_data_info_sptr())); } - scatter_estimate_sptr->fill(0.0); + scatter_estimate_sptr->fill(0.0); // Upsample to 3D //we're currently not doing the tail fitting in this step, but keeping the same scale as determined in 2D - //Note that most of the arguments here are ignored because we fix the scale to 1 - shared_ptr normalisation_factors_3d_sptr = - this->get_normalisation_object_sptr(this->multiplicative_binnorm_sptr); + //Note that most of the arguments here are ignored because we fix the scale to 1 + shared_ptr normalisation_factors_3d_sptr = + this->get_normalisation_object_sptr(this->multiplicative_binnorm_sptr); upsample_and_fit_scatter_estimate(*scatter_estimate_sptr, *this->input_projdata_sptr, - *temp_projdata, + temp_projdata, *normalisation_factors_3d_sptr, *this->input_projdata_sptr, 1.0f, 1.0f, 1, spline_type, - false); - } - else - { - scatter_estimate_sptr = scaled_est_projdata_sptr; - } - - if(!this->output_additive_estimate_prefix.empty()) - { - info("ScatterEstimation: constructing additive sinogram"); - // Now save the full background term. - std::stringstream convert; - convert << this->output_additive_estimate_prefix << "_" << - i_scat_iter; - std::string output_additive_filename = convert.str(); - - shared_ptr temp_additive_projdata( - new ProjDataInterfile(this->input_projdata_sptr->get_exam_info_sptr(), - this->input_projdata_sptr->get_proj_data_info_sptr() , - output_additive_filename, - std::ios::in | std::ios::out | std::ios::trunc)); - - temp_additive_projdata->fill(*scatter_estimate_sptr); - if (!is_null_ptr(this->back_projdata_sptr)) - { - add_proj_data(*temp_additive_projdata, *this->back_projdata_sptr); - } - - this->multiplicative_binnorm_sptr->apply(*temp_additive_projdata); - } + false, false); + } + else + { + scatter_estimate_sptr = scaled_est_projdata_sptr; + } + + if(!this->output_additive_estimate_prefix.empty()) + { + info("ScatterEstimation: constructing additive sinogram"); + // Now save the full background term. + std::stringstream convert; + convert << this->output_additive_estimate_prefix << "_" << + i_scat_iter; + std::string output_additive_filename = convert.str(); + + shared_ptr temp_additive_projdata( + new ProjDataInterfile(this->input_projdata_sptr->get_exam_info_sptr(), + this->input_projdata_sptr->get_proj_data_info_sptr() , + output_additive_filename, + std::ios::in | std::ios::out | std::ios::trunc)); + + temp_additive_projdata->fill(*scatter_estimate_sptr); + if (!is_null_ptr(this->back_projdata_sptr)) + { + add_proj_data(*temp_additive_projdata, *this->back_projdata_sptr); + } + + // Apply * + this->multiplicative_binnorm_sptr->apply(*temp_additive_projdata); + } } // In the additive put the scaled scatter estimate @@ -1090,7 +1114,14 @@ process_data() else { // TODO restructure code to move additive_projdata code from above - error("ScatterEstimation: You should not be here. This is not 2D."); + warning("ScatterEstimation: Great this is 3D!"); + this->add_projdata_sptr->fill(*scaled_est_projdata_sptr); + + if (!is_null_ptr(this->back_projdata_sptr)) + { + add_proj_data(*add_projdata_sptr, *this->back_projdata_sptr); + } + this->multiplicative_binnorm_sptr->apply(*add_projdata_sptr); } current_activity_image_sptr->fill(1.f); diff --git a/src/scatter_buildblock/upsample_and_fit_scatter_estimate.cxx b/src/scatter_buildblock/upsample_and_fit_scatter_estimate.cxx index e002d45586..2295a6833e 100644 --- a/src/scatter_buildblock/upsample_and_fit_scatter_estimate.cxx +++ b/src/scatter_buildblock/upsample_and_fit_scatter_estimate.cxx @@ -3,8 +3,8 @@ Copyright (C) 2014, 2020 University College London This file is part of STIR. - SPDX-License-Identifier: Apache-2.0 - + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for details */ /*! @@ -26,7 +26,7 @@ #include "stir/recon_buildblock/BinNormalisation.h" #include "stir/interpolate_projdata.h" #include "stir/utilities.h" -#include "stir/IndexRange2D.h" +#include "stir/IndexRange2D.h" #include "stir/stream.h" #include "stir/Succeeded.h" #include "stir/thresholding.h" @@ -39,124 +39,142 @@ START_NAMESPACE_STIR -void +void ScatterEstimation:: upsample_and_fit_scatter_estimate(ProjData& scaled_scatter_proj_data, const ProjData& emission_proj_data, - const ProjData& scatter_proj_data, + std::shared_ptr< ProjData> scatter_proj_data, BinNormalisation& scatter_normalisation, const ProjData& weights_proj_data, const float min_scale_factor, const float max_scale_factor, const unsigned half_filter_width, BSpline::BSplineType spline_type, - const bool remove_interleaving) + const bool remove_interleaving, + const bool do_3D) { - shared_ptr - interpolated_direct_scatter_proj_data_info_sptr(emission_proj_data.get_proj_data_info_sptr()->clone()); + shared_ptr + interpolated_scatter_proj_data_info_sptr(emission_proj_data.get_proj_data_info_sptr()->clone()); - info("upsample_and_fit_scatter_estimate: Interpolating scatter estimate to size of emission data"); - ProjDataInMemory interpolated_direct_scatter(emission_proj_data.get_exam_info_sptr(), - interpolated_direct_scatter_proj_data_info_sptr); - - bool actual_remove_interleaving = remove_interleaving; + info("upsample_and_fit_scatter_estimate: Interpolating scatter estimate to size of emission data"); + shared_ptr interpolated_scatter( new ProjDataInMemory(emission_proj_data.get_exam_info_sptr(), + interpolated_scatter_proj_data_info_sptr)); - if (remove_interleaving && emission_proj_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_scanner_geometry()!="Cylindrical") - { - warning("upsample_and_fit_scatter_estimate: forcing remove_interleaving to false as non-cylindrical projdata"); - actual_remove_interleaving = false; - } - if (emission_proj_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_scanner_geometry()=="Cylindrical") - interpolated_direct_scatter_proj_data_info_sptr->reduce_segment_range(0,0); + bool actual_remove_interleaving = remove_interleaving; - interpolate_projdata(interpolated_direct_scatter, scatter_proj_data, spline_type, actual_remove_interleaving); + if (remove_interleaving && emission_proj_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_scanner_geometry()!="Cylindrical") + { + warning("upsample_and_fit_scatter_estimate: forcing remove_interleaving to false as non-cylindrical projdata, or trying 3D scatter estimation"); + actual_remove_interleaving = false; + } + if (emission_proj_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_scanner_geometry()=="Cylindrical" && !do_3D) + { + interpolated_scatter_proj_data_info_sptr->reduce_segment_range(0,0); + interpolate_projdata(*interpolated_scatter, scatter_proj_data, + spline_type, actual_remove_interleaving); + } + else + interpolate_projdata_3d(interpolated_scatter, scatter_proj_data, + spline_type); - const TimeFrameDefinitions& time_frame_defs = - emission_proj_data.get_exam_info_sptr()->time_frame_definitions; + const TimeFrameDefinitions& time_frame_defs = + emission_proj_data.get_exam_info_sptr()->time_frame_definitions; - if (min_scale_factor != 1 || max_scale_factor != 1 || !scatter_normalisation.is_trivial()) + if (min_scale_factor != 1 || max_scale_factor != 1 || !scatter_normalisation.is_trivial()) { - ProjDataInMemory interpolated_scatter(emission_proj_data.get_exam_info_sptr(), - emission_proj_data.get_proj_data_info_sptr()->create_shared_clone()); - if (emission_proj_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_scanner_geometry()!="Cylindrical") - { - interpolated_scatter.fill_from(interpolated_direct_scatter.begin()); - } - else - inverse_SSRB(interpolated_scatter, interpolated_direct_scatter); - - scatter_normalisation.set_up(emission_proj_data.get_exam_info_sptr(), emission_proj_data.get_proj_data_info_sptr()->create_shared_clone()); - scatter_normalisation.undo(interpolated_scatter); - Array<2,float> scale_factors; - - if (min_scale_factor == max_scale_factor) - { - if (min_scale_factor == 1.F) + shared_ptr _interpolated_scatter; + if(!do_3D) + { + _interpolated_scatter.reset(new ProjDataInMemory(emission_proj_data.get_exam_info_sptr(), + emission_proj_data.get_proj_data_info_sptr()->create_shared_clone())); + if (emission_proj_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_scanner_geometry()!="Cylindrical") + { + _interpolated_scatter->fill_from(interpolated_scatter->begin()); + } + else + { + inverse_SSRB(*_interpolated_scatter, *interpolated_scatter); + } + } + else + _interpolated_scatter = interpolated_scatter; + scatter_normalisation.set_up(emission_proj_data.get_exam_info_sptr(), emission_proj_data.get_proj_data_info_sptr()->create_shared_clone()); + // div + scatter_normalisation.undo(*_interpolated_scatter); + Array<2,float> scale_factors; + + if (min_scale_factor == max_scale_factor) + { + if (min_scale_factor == 1.F) { - scaled_scatter_proj_data.fill(interpolated_scatter); - return; // all done + scaled_scatter_proj_data.fill(*_interpolated_scatter); + return; // all done } - const ProjDataInfo& proj_data_info = *emission_proj_data.get_proj_data_info_sptr(); - IndexRange2D sinogram_range(proj_data_info.get_min_segment_num(),proj_data_info.get_max_segment_num(),0,0); - for (int segment_num=proj_data_info.get_min_segment_num(); - segment_num<=proj_data_info.get_max_segment_num(); - ++segment_num) - { - sinogram_range[segment_num].resize( - proj_data_info.get_min_axial_pos_num(segment_num), - proj_data_info.get_max_axial_pos_num(segment_num) ); - } - scale_factors.grow(sinogram_range); - scale_factors.fill(min_scale_factor); - } - else - { - info("upsample_and_fit_scatter_estimate: Finding scale factors by sinogram", 3); - scale_factors = get_scale_factors_per_sinogram( - emission_proj_data, - interpolated_scatter, - weights_proj_data); - - info(boost::format("upsample_and_fit_scatter_estimate: scale factors before thresholding:\n%1%") % - scale_factors, - 2); - - threshold_lower(scale_factors.begin_all(), - scale_factors.end_all(), - min_scale_factor); - threshold_upper(scale_factors.begin_all(), - scale_factors.end_all(), - max_scale_factor); - info(boost::format("upsample_and_fit_scatter_estimate: scale factors after thresholding:\n%1%") % - scale_factors, - 2); - VectorWithOffset kernel(-static_cast(half_filter_width),half_filter_width); - kernel.fill(1.F/(2*half_filter_width+1)); - ArrayFilter1DUsingConvolution lowpass_filter(kernel, BoundaryConditions::constant); - std::for_each(scale_factors.begin(), - scale_factors.end(), - lowpass_filter); - info(boost::format("upsample_and_fit_scatter_estimate: scale factors after filtering:\n%1%") % - scale_factors, - 2); - } - info("upsample_and_fit_scatter_estimate: applying scale factors", 3); - if (scale_sinograms(scaled_scatter_proj_data, - interpolated_scatter, - scale_factors) != Succeeded::yes) + const ProjDataInfo& proj_data_info = *emission_proj_data.get_proj_data_info_sptr(); + IndexRange2D sinogram_range(proj_data_info.get_min_segment_num(),proj_data_info.get_max_segment_num(),0,0); + for (int segment_num=proj_data_info.get_min_segment_num(); + segment_num<=proj_data_info.get_max_segment_num(); + ++segment_num) + { + sinogram_range[segment_num].resize( + proj_data_info.get_min_axial_pos_num(segment_num), + proj_data_info.get_max_axial_pos_num(segment_num) ); + } + scale_factors.grow(sinogram_range); + scale_factors.fill(min_scale_factor); + } + else { - error("upsample_and_fit_scatter_estimate: writing of scaled sinograms failed"); + info("upsample_and_fit_scatter_estimate: Finding scale factors by sinogram", 3); + scale_factors = get_scale_factors_per_sinogram( + emission_proj_data, + *_interpolated_scatter, + weights_proj_data); + + info(boost::format("upsample_and_fit_scatter_estimate: scale factors before thresholding:\n%1%") % + scale_factors, + 2); + + threshold_lower(scale_factors.begin_all(), + scale_factors.end_all(), + min_scale_factor); + threshold_upper(scale_factors.begin_all(), + scale_factors.end_all(), + max_scale_factor); + info(boost::format("upsample_and_fit_scatter_estimate: scale factors after thresholding:\n%1%") % + scale_factors, + 2); + VectorWithOffset kernel(-static_cast(half_filter_width),half_filter_width); + kernel.fill(1.F/(2*half_filter_width+1)); + ArrayFilter1DUsingConvolution lowpass_filter(kernel, BoundaryConditions::constant); + std::for_each(scale_factors.begin(), + scale_factors.end(), + lowpass_filter); + info(boost::format("upsample_and_fit_scatter_estimate: scale factors after filtering:\n%1%") % + scale_factors, + 2); + } + info("upsample_and_fit_scatter_estimate: applying scale factors", 3); + if (scale_sinograms(scaled_scatter_proj_data, + *_interpolated_scatter, + scale_factors) != Succeeded::yes) + { + error("upsample_and_fit_scatter_estimate: writing of scaled sinograms failed"); } } - else // min/max_scale_factor equal to 1 and no norm + else // min/max_scale_factor equal to 1 and no norm -> Go to 3D if not { - if (emission_proj_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_scanner_geometry()!="Cylindrical") - scaled_scatter_proj_data.fill_from(interpolated_direct_scatter.begin()); - else - inverse_SSRB(scaled_scatter_proj_data, interpolated_direct_scatter); + if (emission_proj_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_scanner_geometry()!="Cylindrical") + scaled_scatter_proj_data.fill_from(interpolated_scatter->begin()); + else + { + if(!do_3D) + inverse_SSRB(scaled_scatter_proj_data, *interpolated_scatter); + } } + std::cout << "Getting out of Upsampling" << std::endl; } END_NAMESPACE_STIR diff --git a/src/scatter_utilities/upsample_and_fit_single_scatter.cxx b/src/scatter_utilities/upsample_and_fit_single_scatter.cxx index 48bc6d6bc5..7f5e17cf87 100644 --- a/src/scatter_utilities/upsample_and_fit_single_scatter.cxx +++ b/src/scatter_utilities/upsample_and_fit_single_scatter.cxx @@ -192,7 +192,7 @@ int main(int argc, const char *argv[]) stir::ScatterEstimation:: upsample_and_fit_scatter_estimate(output_proj_data, *data_to_fit_proj_data_sptr, - *data_to_scale_proj_data_sptr, + data_to_scale_proj_data_sptr, *normalisation_sptr, *weights_proj_data_sptr, min_scale_factor,