ProjMatrixByBinUsingInterpolation.h

Go to the documentation of this file.

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

    SPDX-License-Identifier: Apache-2.0

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

#include "stir/RegisteredParsingObject.h"
#include "stir/recon_buildblock/ProjMatrixByBin.h"
#include "stir/ProjDataInfoCylindrical.h"
#include "stir/CartesianCoordinate3D.h"
#include "stir/shared_ptr.h"

START_NAMESPACE_STIR

template <int num_dimensions, typename elemT>
class DiscretisedDensity;
class Bin;
class ProjMatrixByBinUsingInterpolation
    : public RegisteredParsingObject<ProjMatrixByBinUsingInterpolation, ProjMatrixByBin, ProjMatrixByBin>
{
public:
  static const char* const registered_name;

  ProjMatrixByBinUsingInterpolation();


  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;

  ProjMatrixByBinUsingInterpolation* clone() const override;

private:
  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;

  // explicitly list necessary members for image details (should use an Info object instead)
  CartesianCoordinate3D<float> voxel_size;
  CartesianCoordinate3D<float> origin;
  IndexRange<3> densel_range;

  shared_ptr<const ProjDataInfo> proj_data_info_ptr;

  // for Jacobian
  const ProjDataInfoCylindrical& proj_data_info_cyl() const
  {
    return static_cast<const ProjDataInfoCylindrical&>(*proj_data_info_ptr);
  }
  class JacobianForIntBP
  {
  private:
    // store some scanner related data to avoid recomputation
    float R2;
    float ring_spacing2;
    bool arccor;
    // 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)

    float backprojection_normalisation;

    bool use_exact_Jacobian_now;

  public:
    // default constructor needed as now member of projector class (better to make set_up)
    JacobianForIntBP() {}
    explicit JacobianForIntBP(const ProjDataInfoCylindrical* proj_data_info_ptr, bool exact);
    // s in mm here!
    float operator()(const float delta, const float s) const
    {
      float tmp;
      if (use_exact_Jacobian_now)
        tmp = 4 * (R2 - s * s);
      else
        tmp = 4 * R2;
      if (!arccor)
        tmp *= sqrt(tmp);
      return (arccor ? tmp : pow(tmp, 1.5F)) / pow(tmp + ring_spacing2 * delta * delta, 1.5F) * backprojection_normalisation;
    }
  };

  JacobianForIntBP jacobian;
  bool use_piecewise_linear_interpolation_now;
  bool use_exact_Jacobian_now;

  void calculate_proj_matrix_elems_for_one_bin(ProjMatrixElemsForOneBin&) const override;

  void set_defaults() override;
  void initialise_keymap() override;
  bool post_processing() override;

  float get_element(const Bin& bin, const CartesianCoordinate3D<float>& densel_ctr) const;

private:
  void find_tang_ax_pos_diff(float& tang_pos_diff,
                             float& ax_pos_diff,
                             const Bin& bin,
                             const CartesianCoordinate3D<float>& point) const;
};

END_NAMESPACE_STIR

#endif