KOSMAPOSLReconstruction.h

Go to the documentation of this file.

//
//
/*

    Copyright (C) 2018 Commonwealth Scientific and Industrial Research Organisation
    Copyright (C) 2018-2019 University of Leeds
    Copyright (C) 2019 University College of London
    Copyright (C) 2019-2021 National Physical Laboratory

    This file is part of STIR.

    SPDX-License-Identifier: Apache-2.0

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

#  include "stir/recon_buildblock/IterativeReconstruction.h"
#  include "stir/RegisteredParsingObject.h"
#  include "stir/OSMAPOSL/OSMAPOSLReconstruction.h"
#  include "stir/CartesianCoordinate3D.h"

START_NAMESPACE_STIR

template <typename TargetT>
class PoissonLogLikelihoodWithLinearModelForMean;

template <typename TargetT>
class KOSMAPOSLReconstruction
    : public RegisteredParsingObject<KOSMAPOSLReconstruction<TargetT>, Reconstruction<TargetT>, OSMAPOSLReconstruction<TargetT>>
{
private:
  typedef RegisteredParsingObject<KOSMAPOSLReconstruction<TargetT>, Reconstruction<TargetT>, OSMAPOSLReconstruction<TargetT>>

      base_type;

public:
  static const char* const registered_name;

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

  KOSMAPOSLReconstruction& get_parameters() { return *this; }

  const KOSMAPOSLReconstruction& get_parameters() const { return *this; }

  // kernel
  const std::vector<std::string> get_anatomical_image_filenames() const;
  const int get_num_neighbours() const;
  const int get_num_non_zero_feat() const;
  const std::vector<double> get_sigma_m() const;
  const double get_sigma_p() const;
  const double get_sigma_dp() const;
  const double get_sigma_dm() const;
  const bool get_only_2D() const;
  const bool get_hybrid() const;
  const int get_freeze_iterative_kernel_at_subiter_num() const;

  std::vector<shared_ptr<TargetT>> get_anatomical_prior_sptrs();

  void set_anatomical_prior_sptr(shared_ptr<TargetT>, int index);
  void set_anatomical_prior_sptr(shared_ptr<TargetT> arg);
  void set_anatomical_image_filename(const std::string&, const int index);
  void set_anatomical_image_filename(const std::string&);

  void set_num_neighbours(const int);
  void set_num_non_zero_feat(const int);
  void set_sigma_m(const double, const int index);
  void set_sigma_m(const double);
  void set_sigma_p(const double);
  void set_sigma_dp(const double);
  void set_sigma_dm(const double);
  void set_only_2D(const bool);
  void set_hybrid(const bool);
  void set_freeze_iterative_kernel_at_subiter_num(const int);

  void ask_parameters() override;

protected:
  std::string input_filename, kernelised_output_filename_prefix;
  std::string current_kimage_filename;
  std::string sens_filenames;

  std::vector<std::string> anatomical_image_filenames;

  std::vector<shared_ptr<TargetT>> anatomical_prior_sptrs, kmnorm_sptrs;
  shared_ptr<TargetT> kpnorm_sptr;
  // kernel parameters
  int num_neighbours, num_non_zero_feat, num_elem_neighbourhood, num_voxels, dimz, dimy, dimx;
  int freeze_iterative_kernel_at_subiter_num;
  std::vector<double> sigma_m;
  bool only_2D;
  bool hybrid;
  double sigma_p;
  double sigma_dp, sigma_dm;
  BasicCoordinate<3, int> min_ind, max_ind;
  shared_ptr<TargetT> iterative_kernel_image_frozen_sptr;

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

  void compute_kernelised_image(TargetT& kernelised_image_out,
                                const TargetT& image_to_kernelise,
                                const TargetT& current_alpha_estimate);

private:
  friend void do_sensitivity(const char* const par_filename);

  Succeeded set_up(shared_ptr<TargetT> const& target_image_sptr) override;

  void update_estimate(TargetT& current_image_estimate) override;

  bool still_updating_iterative_kernel();

  std::vector<double> anatomical_sd;
  mutable Array<3, float> distance;
  void calculate_norm_matrix(TargetT& normp, const int dimf_row, const int dimf_col, const TargetT& emission);

  void calculate_norm_const_matrix(std::vector<shared_ptr<TargetT>>& normm, const int dimf_row, const int dimf_col);

  void estimate_stand_dev_for_anatomical_image(std::vector<double>& SD);
#  if 0

//    void full_compute_kernelised_image(TargetT& kernelised_image_out,
//                                     const TargetT& image_to_kernelise,
//                                     const TargetT& current_alpha_estimate);

//    void compact_compute_kernelised_image(TargetT& kernelised_image_out,
//                                        const TargetT& image_to_kernelise,
//                                        const TargetT& current_alpha_estimate);
#  endif

  double calc_emission_kernel(const double current_alpha_estimate_zyx,
                              const double current_alpha_estimate_zyx_dr,
                              const double distance_dzdydx,
                              const bool use_compact_implementation,
                              const int l,
                              const int m);

  double calc_anatomical_kernel(const double anatomical_prior_zyx,
                                const double anatomical_prior_zyx_dr,
                                const double distance_dzdydx,
                                const bool use_compact_implementation,
                                const int l,
                                const int m,
                                const int index);

  double calc_kernel_from_precalculated(const double precalculated_norm_zxy,
                                        const double sq_sigma_int,
                                        const double sq_sigma_dist,
                                        const double sq_distance_dzdydx,
                                        const double precalc_denom);

  double calc_kernel_compact(const double prior_image_zyx_diff,
                             const double sq_sigma_int,
                             const double sq_sigma_dist,
                             const double sq_distance_dzdydx,
                             const double precalc_denom);
};

END_NAMESPACE_STIR

#endif

// __KOSMAPOSLReconstruction_h__