ScatterEstimation.h

Go to the documentation of this file.

#ifndef __stir_scatter_ScatterEstimation_H__
#define __stir_scatter_ScatterEstimation_H__

/*
    Copyright (C) 2018 - 2019 University of Hull
    Copyright (C) 2016,2020 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
*/
#include "stir/DiscretisedDensity.h"
#include "stir/ProjData.h"
#include "stir/ParsingObject.h"
#include "stir/numerics/BSplines.h"
#include "stir/CartesianCoordinate3D.h"
#include "stir/recon_buildblock/Reconstruction.h"

#include "stir/scatter/ScatterSimulation.h"

#include "stir/recon_buildblock/IterativeReconstruction.h"
#include "stir/recon_buildblock/GeneralisedObjectiveFunction.h"
#include "stir/recon_buildblock/AnalyticReconstruction.h"

#include "stir/stir_math.h"
#include "stir/FilePath.h"
#include "stir/recon_buildblock/ForwardProjectorByBin.h"

START_NAMESPACE_STIR

template <class TargetT>
class PostFiltering;
class BinNormalisation;

struct MaskingParameters
{
  float min_threshold;
  std::string filter_filename;
  shared_ptr<PostFiltering<DiscretisedDensity<3, float>>> filter_sptr;
};

class ScatterEstimation : public ParsingObject
{
public:

  static void upsample_and_fit_scatter_estimate(ProjData& scaled_scatter_proj_data,
                                                const ProjData& emission_proj_data,
                                                const 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 = BSpline::BSplineType::linear,
                                                const bool remove_interleaving = true);

  ScatterEstimation();
  explicit ScatterEstimation(const std::string& parameter_filename);

  virtual Succeeded process_data();

  shared_ptr<ProjData> get_output() const;

  shared_ptr<ProjData> make_2D_projdata_sptr(const shared_ptr<ProjData> in_3d_sptr);
  shared_ptr<ProjData>
  make_2D_projdata_sptr(const shared_ptr<ProjData> in_3d_sptr, string template_filename, const bool do_normalisation = false);

  virtual Succeeded set_up();

  // Set functions
  inline void set_input_proj_data_sptr(const shared_ptr<ProjData>);

#if STIR_VERSION < 050000
  void set_input_data(const shared_ptr<ProjData>& data);
#else
  void set_input_data(const shared_ptr<ExamData>& data);
#endif
  shared_ptr<const ProjData> get_input_data() const;

  inline void set_reconstruction_method_sptr(const shared_ptr<Reconstruction<DiscretisedDensity<3, float>>>);
  inline void set_attenuation_image_sptr(const shared_ptr<const DiscretisedDensity<3, float>>);
  void set_attenuation_correction_proj_data_sptr(const shared_ptr<ProjData>);
  void set_normalisation_sptr(const shared_ptr<BinNormalisation>);
  inline void set_background_proj_data_sptr(const shared_ptr<ProjData>);
  inline void set_initial_activity_image_sptr(const shared_ptr<const DiscretisedDensity<3, float>>);

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

  inline void set_mask_proj_data_sptr(const shared_ptr<ProjData> arg);

  void set_recompute_mask_image(bool arg);
  void set_recompute_mask_projdata(bool arg);

  void set_forward_projector_for_mask_sptr(const shared_ptr<ForwardProjectorByBin> projector_sptr);
  shared_ptr<ForwardProjectorByBin> get_forward_projector_for_mask_sptr() const;

  inline void set_scatter_simulation_method_sptr(const shared_ptr<ScatterSimulation>);

  inline void set_num_iterations(int);

  inline unsigned int get_half_filter_width() const;
  inline void set_half_filter_width(unsigned int);

  inline void
  set_downsample_scanner(bool downsample_scanner, int downsampled_number_of_rings = -1, int downsampled_detectors_per_ring = -1);

  void set_output_scatter_estimate_prefix(const std::string&);
  void set_export_scatter_estimates_of_each_iteration(bool);

  void set_max_scale_value(float value);
  void set_min_scale_value(float value);
  void set_mask_projdata_filename(std::string name);
  void set_mask_image_filename(std::string name);
  void set_output_additive_estimate_prefix(std::string name);
  void set_run_debug_mode(bool debug);
  void set_restart_reconstruction_every_scatter_iteration(bool setting);
  bool get_restart_reconstruction_every_scatter_iteration() const;

  // inline void set_zoom_xy(float);
  // inline void set_zoom_z(float);

  // Get functions

  int get_iterations_num() const;

  int get_num_iterations() const;

  shared_ptr<const DiscretisedDensity<3, float>> get_estimated_activity_image_sptr() const;

  virtual bool already_setup() const;

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

  bool recompute_mask_image;
  bool recompute_mask_projdata;
  bool recompute_atten_projdata;
  bool restart_reconstruction_every_scatter_iteration;

  shared_ptr<Reconstruction<DiscretisedDensity<3, float>>> reconstruction_template_sptr;
  shared_ptr<DiscretisedDensity<3, float>> current_activity_image_sptr;
  shared_ptr<const DiscretisedDensity<3, float>> atten_image_sptr;
  shared_ptr<BinNormalisation> norm_3d_sptr;
  shared_ptr<ProjData> mask_projdata_sptr;
  shared_ptr<ProjData> input_projdata_sptr;
  shared_ptr<ProjData> input_projdata_2d_sptr;
  shared_ptr<ProjData> add_projdata_2d_sptr;
  shared_ptr<ProjData> data_to_fit_projdata_sptr;

  shared_ptr<ProjData> add_projdata_sptr;
  shared_ptr<ProjData> back_projdata_2d_sptr;
  shared_ptr<ProjData> back_projdata_sptr;
  std::string mask_image_filename;
  std::string mask_projdata_filename;
  std::string back_projdata_filename;

  std::string tail_mask_par_filename;
  std::string input_projdata_filename;
  std::string atten_image_filename;
  std::string recon_template_par_filename;
  std::string atten_coeff_filename;
  shared_ptr<ForwardProjectorByBin> forward_projector_for_mask_sptr;


  MaskingParameters masking_parameters;
  int num_scatter_iterations;
  std::string output_scatter_estimate_prefix;

  std::string output_additive_estimate_prefix;

  bool _already_setup;

private:
  shared_ptr<BinNormalisation> atten_norm_3d_sptr;


  shared_ptr<BinNormalisation> multiplicative_binnorm_2d_sptr;


  shared_ptr<BinNormalisation> multiplicative_binnorm_sptr;

  shared_ptr<ProjData> scatter_estimate_sptr;

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

  Succeeded set_up_iterative(shared_ptr<IterativeReconstruction<DiscretisedDensity<3, float>>> arg);

  Succeeded set_up_analytic();

  Succeeded project_mask_image();


  void reconstruct_iterative(int scat_iter);


  void reconstruct_analytic(int scat_iter);

  static void apply_mask_in_place(DiscretisedDensity<3, float>&, const MaskingParameters&);

  void add_proj_data(ProjData&, const ProjData&);

  void subtract_proj_data(ProjData&, const ProjData&);

  void apply_to_proj_data(ProjData&, const pow_times_add&);

  void create_multiplicative_binnorm_sptr();

  shared_ptr<BinNormalisation> get_normalisation_object_sptr(const shared_ptr<BinNormalisation>& combined_norm_sptr) const;

  shared_ptr<ProjData> get_attenuation_correction_factors_sptr(const shared_ptr<BinNormalisation>& combined_norm_sptr) const;

  shared_ptr<ProjData> create_new_proj_data(const std::string& filename,
                                            const shared_ptr<const ExamInfo> exam_info_sptr,
                                            const shared_ptr<const ProjDataInfo> proj_data_info_sptr) const;

  bool do_average_at_2;
  bool remove_interleaving;
  bool export_scatter_estimates_of_each_iteration;
  bool run_in_2d_projdata;
  bool override_density_image;
  bool override_scanner_template;
  bool run_debug_mode;
  std::string scatter_sim_par_filename;
  shared_ptr<ScatterSimulation> scatter_simulation_sptr;
  FilePath extras_path;

  float max_scale_value;
  float min_scale_value;

  bool downsample_scanner_bool;
  int downsampled_number_of_rings;
  int downsampled_detectors_per_ring;

  unsigned int half_filter_width;

  bool iterative_method;
};

END_NAMESPACE_STIR
#include "stir/scatter/ScatterEstimation.inl"
#endif