ScatterSimulation.h

Go to the documentation of this file.

/*
    Copyright (C) 2018 - 2019 University of Hull
    Copyright (C) 2004 - 2009 Hammersmith Imanet Ltd
    Copyright (C) 2013 - 2016, 2019, 2020, 2022 University College London
    Copyright (C) 2022, National Physical Laboratory
    This file is part of STIR.

    SPDX-License-Identifier: Apache-2.0

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

#include "stir/ParsingObject.h"
#include "stir/RegisteredObject.h"
#include "stir/ProjData.h"
#include "stir/VoxelsOnCartesianGrid.h"
#include "stir/ProjDataInfoBlocksOnCylindricalNoArcCorr.h"
#include "stir/ProjDataInfoCylindricalNoArcCorr.h"
#include "stir/ProjDataInfoGenericNoArcCorr.h"

START_NAMESPACE_STIR

class ScatterSimulation : public RegisteredObject<ScatterSimulation>
{
public:
  ScatterSimulation();

  ~ScatterSimulation() override;

  virtual Succeeded process_data();
  virtual std::string method_info() const = 0;
  virtual void ask_parameters();

  inline bool has_template_proj_data_info() const { return !stir::is_null_ptr(proj_data_info_sptr); }
  inline bool has_exam_info() const { return !stir::is_null_ptr(template_exam_info_sptr); }


  shared_ptr<ProjData> get_output_proj_data_sptr() const;

  inline int get_num_scatter_points() const { return static_cast<int>(this->scatt_points_vector.size()); }
  shared_ptr<const ProjDataInfo> get_template_proj_data_info_sptr() const;
  shared_ptr<const ExamInfo> get_exam_info_sptr() const;

  const DiscretisedDensity<3, float>& get_activity_image() const;
  const DiscretisedDensity<3, float>& get_attenuation_image() const;
  const DiscretisedDensity<3, float>& get_attenuation_image_for_scatter_points() const;
  shared_ptr<const DiscretisedDensity<3, float>> get_density_image_for_scatter_points_sptr() const;



  void set_template_proj_data_info(const std::string&);

  void set_template_proj_data_info(const ProjDataInfo&);

  void set_activity_image_sptr(const shared_ptr<const DiscretisedDensity<3, float>>);

  void set_activity_image(const std::string& filename);
  void set_exam_info(const ExamInfo&);
  void set_exam_info_sptr(const shared_ptr<const ExamInfo>);

  void set_output_proj_data_sptr(shared_ptr<ProjData>);

  void set_density_image_sptr(const shared_ptr<const DiscretisedDensity<3, float>>);

  void set_density_image(const std::string&);
  void set_output_proj_data(const std::string&);

  void set_output_proj_data_sptr(const shared_ptr<const ExamInfo>, const shared_ptr<const ProjDataInfo>, const std::string&);

  void set_density_image_for_scatter_points_sptr(shared_ptr<const DiscretisedDensity<3, float>>);

  void set_image_downsample_factors(float factor_xy = 1.f, float factor_z = 1.f, int _size_zoom_xy = -1, int _size_zoom_z = -1);
  void set_density_image_for_scatter_points(const std::string&);
  void set_attenuation_threshold(const float);

  void set_randomly_place_scatter_points(const bool);

  void set_cache_enabled(const bool);


  void downsample_density_image_for_scatter_points(float _zoom_xy, float _zoom_z, int _size_xy = -1, int _size_z = -1);


  void set_downsample_scanner_bool(const bool arg);
  bool get_downsample_scanner_bool() const;


  int get_num_downsample_scanner_rings() const;
  void set_num_downsample_scanner_rings(const int arg);


  int get_num_downsample_scanner_dets() const;
  void set_num_downsample_scanner_dets(const int arg);


  Succeeded downsample_scanner(int new_num_rings = -1, int new_num_dets = -1);

  Succeeded downsample_images_to_scanner_size();


  float detection_efficiency(const float energy) const;


  static inline float dif_Compton_cross_section(const float cos_theta, float energy);

  static inline float total_Compton_cross_section(float energy);

  static inline float photon_energy_after_Compton_scatter(const float cos_theta, const float energy);

  static inline float photon_energy_after_Compton_scatter_511keV(const float cos_theta);

  static inline float total_Compton_cross_section_relative_to_511keV(const float energy);

  virtual Succeeded set_up();

  virtual void write_log(const double simulation_time, const float total_scatter);

  void set_use_cache(const bool);
  bool get_use_cache() const;

protected:

  virtual double process_data_for_view_segment_num(const ViewSegmentNumbers& vs_num);

  float compute_emis_to_det_points_solid_angle_factor(const CartesianCoordinate3D<float>& emis_point,
                                                      const CartesianCoordinate3D<float>& detector_coord);

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

  enum image_type
  {
    act_image_type,
    att_image_type
  };
  struct ScatterPoint
  {
    CartesianCoordinate3D<float> coord;
    float mu_value;
  };

  std::vector<ScatterPoint> scatt_points_vector;

  float scatter_volume;


  void sample_scatter_points();


  virtual void remove_cache_for_integrals_over_attenuation();


  virtual void remove_cache_for_integrals_over_activity();

  static float max_cos_angle(const float low, const float approx, const float resolution_at_511keV);

  static float energy_lower_limit(const float low, const float approx, const float resolution_at_511keV);

  virtual void find_detectors(unsigned& det_num_A, unsigned& det_num_B, const Bin& bin) const;

  unsigned find_in_detection_points_vector(const CartesianCoordinate3D<float>& coord) const;

  CartesianCoordinate3D<float> shift_detector_coordinates_to_origin;

  double detection_efficiency_no_scatter(const unsigned det_num_A, const unsigned det_num_B) const;

  // next needs to be mutable because find_in_detection_points_vector is const
  mutable std::vector<CartesianCoordinate3D<float>> detection_points_vector;


  virtual double scatter_estimate(const Bin& bin) = 0;


  static float integral_between_2_points(const DiscretisedDensity<3, float>& density,
                                         const CartesianCoordinate3D<float>& point1,
                                         const CartesianCoordinate3D<float>& point2);

  float exp_integral_over_attenuation_image_between_scattpoint_det(const CartesianCoordinate3D<float>& scatter_point,
                                                                   const CartesianCoordinate3D<float>& detector_coord);

  float integral_over_activity_image_between_scattpoint_det(const CartesianCoordinate3D<float>& scatter_point,
                                                            const CartesianCoordinate3D<float>& detector_coord);

  float cached_integral_over_activity_image_between_scattpoint_det(const unsigned scatter_point_num, const unsigned det_num);

  float cached_exp_integral_over_attenuation_image_between_scattpoint_det(const unsigned scatter_point_num,
                                                                          const unsigned det_num);

  std::string template_proj_data_filename;

  shared_ptr<ProjDataInfo> proj_data_info_sptr;
  shared_ptr<ExamInfo> template_exam_info_sptr;

  std::string density_image_filename;

  std::string density_image_for_scatter_points_filename;

  std::string density_image_for_scatter_points_output_filename;

  shared_ptr<const DiscretisedDensity<3, float>> density_image_sptr;

  shared_ptr<const DiscretisedDensity<3, float>> activity_image_sptr;


  void initialise_cache_for_scattpoint_det_integrals_over_attenuation();

  void initialise_cache_for_scattpoint_det_integrals_over_activity();

  std::string output_proj_data_filename;
  shared_ptr<ProjData> output_proj_data_sptr;

  float attenuation_threshold;

  bool randomly_place_scatter_points;

  bool use_cache;
  std::string activity_image_filename;
  float zoom_xy;
  float zoom_z;
  int zoom_size_xy;
  int zoom_size_z;
  int downsample_scanner_rings;
  int downsample_scanner_dets;

  bool downsample_scanner_bool;
  bool _already_set_up;

private:
  int total_detectors;

  Array<2, float> cached_activity_integral_scattpoint_det;
  Array<2, float> cached_attenuation_integral_scattpoint_det;
  shared_ptr<DiscretisedDensity<3, float>> density_image_for_scatter_points_sptr;

  // numbers that we don't want to recompute all the time
  mutable float detector_efficiency_no_scatter;


  void check_z_to_middle_consistent(const DiscretisedDensity<3, float>& _image, const std::string& name) const;
};

END_NAMESPACE_STIR

#include "stir/scatter/ScatterSimulation.inl"

#endif