BackProjectorByBinUsingInterpolation.h

Go to the documentation of this file.

//
//
/*
    Copyright (C) 2000 PARAPET partners
    Copyright (C) 2000- 2008, Hammersmith Imanet Ltd
    This file is part of STIR.

    SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license

    See STIR/LICENSE.txt for details
*/
#ifndef __BackProjectorByBinUsingInterpolation_h_
#define __BackProjectorByBinUsingInterpolation_h_

#include "stir/recon_buildblock/BackProjectorByBin.h"
#include "stir/RegisteredParsingObject.h"
#include "stir/shared_ptr.h"
#include "stir/ArrayFwd.h"

START_NAMESPACE_STIR

template <typename elemT>
class Viewgram;
template <typename elemT>
class RelatedViewgrams;
template <typename elemT>
class VoxelsOnCartesianGrid;
class ProjDataInfo;
class ProjDataInfoCylindricalArcCorr;
class DataSymmetriesForBins_PET_CartesianGrid;
class JacobianForIntBP
{
private:
  // store some scanner related data to avoid recomputation
  const float R2;
  const float dxy2;
  const float ring_spacing2;
  // total normalisation of backprojection, 3 factors:
  //  (_Pi/scanner.num_views) for discretisation of integral over phi
  // scanner.ring_spacing for discretisation of integral over delta
  // normalisation of projection space integral: 1/(2 Pi)

  const float backprojection_normalisation;

  const bool use_exact_Jacobian_now;

public:
  explicit JacobianForIntBP(const shared_ptr<const ProjDataInfoCylindricalArcCorr> proj_data_info_sptr, bool exact);

  float operator()(const float delta, const float s) const
  {
    float tmp;
    if (use_exact_Jacobian_now)
      tmp = 4 * (R2 - dxy2 * s * s);
    else
      tmp = 4 * R2;
    return tmp / pow(tmp + ring_spacing2 * delta * delta, 1.5F) * backprojection_normalisation;
  }
};

class BackProjectorByBinUsingInterpolation
    : public RegisteredParsingObject<BackProjectorByBinUsingInterpolation, BackProjectorByBin>

{
public:
  static const char* const registered_name;

  explicit BackProjectorByBinUsingInterpolation(const bool use_piecewise_linear_interpolation = true,
                                                const bool use_exact_Jacobian = true);


  BackProjectorByBinUsingInterpolation(shared_ptr<const ProjDataInfo> const&,
                                       shared_ptr<const DiscretisedDensity<3, float>> const& image_info_ptr,
                                       const bool use_piecewise_linear_interpolation = true,
                                       const bool use_exact_Jacobian = true);


  void set_up(const shared_ptr<const ProjDataInfo>& proj_data_info_ptr,
              const shared_ptr<const DiscretisedDensity<3, float>>& density_info_ptr // TODO should be Info only
              ) override;

  const DataSymmetriesForViewSegmentNumbers* get_symmetries_used() const override;
  void use_exact_Jacobian(const bool use_exact_Jacobian);

  void use_piecewise_linear_interpolation(const bool use_piecewise_linear_interpolation);

  BackProjectorByBinUsingInterpolation* clone() const override;

private:
  // KT 20/06/2001 changed type to enable use of more methods
  shared_ptr<DataSymmetriesForBins_PET_CartesianGrid> symmetries_ptr;
  // const DataSymmetriesForViewSegmentNumbers * symmetries_ptr;

  bool use_piecewise_linear_interpolation_now;

  bool use_exact_Jacobian_now;


  bool do_symmetry_90degrees_min_phi;
  bool do_symmetry_180degrees_min_phi;
  bool do_symmetry_swap_segment;
  bool do_symmetry_swap_s;
  bool do_symmetry_shift_z;

#if 0
  // not used yet
struct ProjDataForIntBP
{ 
  float view__pos_s;
  float view__pos_sp1;
  float view__neg_s; 
  float view__neg_sp1;
  float min90__pos_s; 
  float min90__pos_sp1;
  float min90__neg_s; 
  float min90__neg_sp1;
  float plus90__pos_s; 
  float plus90__pos_sp1;
  float plus90__neg_s; 
  float plus90__neg_sp1;
  float min180__pos_s; 
  float min180__pos_sp1;
  float min180__neg_s; 
  float min180__neg_sp1;
};
#endif

  void actual_back_project(DiscretisedDensity<3, float>&,
                           const RelatedViewgrams<float>&,
                           const int min_axial_pos_num,
                           const int max_axial_pos_num,
                           const int min_tangential_pos_num,
                           const int max_tangential_pos_num) override;

  void actual_back_project(DiscretisedDensity<3, float>&, const Bin&);

  virtual void back_project_all_symmetries(VoxelsOnCartesianGrid<float>& image,
                                           const Viewgram<float>& pos_view,
                                           const Viewgram<float>& neg_view,
                                           const Viewgram<float>& pos_plus90,
                                           const Viewgram<float>& neg_plus90,
                                           const Viewgram<float>& pos_min180,
                                           const Viewgram<float>& neg_min180,
                                           const Viewgram<float>& pos_min90,
                                           const Viewgram<float>& neg_min90,
                                           const int min_axial_pos_num,
                                           const int max_axial_pos_num,
                                           const int min_tangential_pos_num,
                                           const int max_tangential_pos_num);

  /*
    This function projects 4 viewgrams related by symmetry.
    It will be used for view=0 or 45 degrees
    (or all others if the above version is not implemented in
    the derived class)
    Here 0<=view < num_views/2 (= 90 degrees)
    */

  virtual void back_project_view_plus_90_and_delta(VoxelsOnCartesianGrid<float>& image,
                                                   const Viewgram<float>& pos_view,
                                                   const Viewgram<float>& neg_view,
                                                   const Viewgram<float>& pos_plus90,
                                                   const Viewgram<float>& neg_plus90,
                                                   const int min_axial_pos_num,
                                                   const int max_axial_pos_num,
                                                   const int min_tangential_pos_num,
                                                   const int max_tangential_pos_num);
  /*
  void back_project_2D_view_plus_90(const PETSinogram<float> &sino, PETPlane &image, int view,
                               const int min_bin_num, const intmax_tangential_pos_num);
  void back_project_2D_all_symmetries(const PETSinogram<float> &sino, PETPlane &image, int view,
                                    const int min_bin_num, const intmax_tangential_pos_num);
*/

  /*

    These functions use a 3D version of Cho's algorithm for backprojecting incrementally.
    See M. Egger's thesis for details.
    In addition to the symmetries mentioned above, they also use s,-s symmetry
    (while being careful when s=0 to avoid self-symmetric cases)
    */

  static void piecewise_linear_interpolation_backproj3D_Cho_view_viewplus90(
      Array<4, float> const& Projptr,
      VoxelsOnCartesianGrid<float>& image,
      const shared_ptr<const ProjDataInfoCylindricalArcCorr> proj_data_info_sptr,
      float delta,
      const double cphi,
      const double sphi,
      int s,
      int ax_pos0,
      const int num_planes_per_axial_pos,
      const float axial_pos_to_z_offset);

  static void piecewise_linear_interpolation_backproj3D_Cho_view_viewplus90_180minview_90minview(
      Array<4, float> const& Projptr,
      VoxelsOnCartesianGrid<float>& image,
      const shared_ptr<const ProjDataInfoCylindricalArcCorr> proj_data_info_sptr,
      float delta,
      const double cphi,
      const double sphi,
      int s,
      int ax_pos0,
      const int num_planes_per_axial_pos,
      const float axial_pos_to_z_offset);

  static void
  linear_interpolation_backproj3D_Cho_view_viewplus90(Array<4, float> const& Projptr,
                                                      VoxelsOnCartesianGrid<float>& image,
                                                      const shared_ptr<const ProjDataInfoCylindricalArcCorr> proj_data_info_sptr,
                                                      float delta,
                                                      const double cphi,
                                                      const double sphi,
                                                      int s,
                                                      int ax_pos0,
                                                      const int num_planes_per_axial_pos,
                                                      const float axial_pos_to_z_offset);

  static void linear_interpolation_backproj3D_Cho_view_viewplus90_180minview_90minview(
      Array<4, float> const& Projptr,
      VoxelsOnCartesianGrid<float>& image,
      const shared_ptr<const ProjDataInfoCylindricalArcCorr> proj_data_info_sptr,
      float delta,
      const double cphi,
      const double sphi,
      int s,
      int ax_pos0,
      const int num_planes_per_axial_pos,
      const float axial_pos_to_z_offset);

  /*
static void   backproj2D_Cho_view_viewplus90( PETPlane & image,
                                     const ProjDataForIntBP &projs,
                                    const double cphi, const double sphi, int s);

  static void   backproj2D_Cho_view_viewplus90_180minview_90minview( PETPlane & image,
                                     const ProjDataForIntBP &projs,
                                    const double cphi, const double sphi, int s);

*/
  void set_defaults() override;
  void initialise_keymap() override;
};

END_NAMESPACE_STIR

#endif // __BackProjectorByBinUsingInterpolation_h_