ParametricQuadraticPrior.h

Go to the documentation of this file.

//
//
/*
    Copyright (C) 2006- 2009, 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_ParametricQuadraticPrior_H__
#define __stir_recon_buildblock_ParametricQuadraticPrior_H__

#include "stir/RegisteredParsingObject.h"
#include "stir/recon_buildblock/QuadraticPrior.h"
#include "stir/recon_buildblock/PriorWithParabolicSurrogate.h"
#include "stir/Array.h"
#include "stir/modelling/ParametricDiscretisedDensity.h"
#include "stir/modelling/KineticParameters.h"
#include "stir/VectorWithOffset.h"
#include "stir/DiscretisedDensity.h"
#include "stir/shared_ptr.h"
#include <string>

START_NAMESPACE_STIR

template <typename TargetT>
class ParametricQuadraticPrior : public RegisteredParsingObject<ParametricQuadraticPrior<TargetT>,
                                                                GeneralisedPrior<TargetT>,
                                                                PriorWithParabolicSurrogate<TargetT>>
{
private:
  typedef RegisteredParsingObject<ParametricQuadraticPrior<TargetT>,
                                  GeneralisedPrior<TargetT>,
                                  PriorWithParabolicSurrogate<TargetT>>
      base_type;
  shared_ptr<TargetT> kappa_ptr;
  VectorWithOffset<QuadraticPrior<float>> _single_quadratic_priors;

public:
  static const char* const registered_name;

  ParametricQuadraticPrior();

  ParametricQuadraticPrior(const bool only_2D, float penalization_factor);

  virtual bool parabolic_surrogate_curvature_depends_on_argument() const { return false; }

  double compute_value(const TargetT& current_image_estimate);

  void compute_gradient(TargetT& prior_gradient, const TargetT& current_image_estimate);


  void parabolic_surrogate_curvature(TargetT& parabolic_surrogate_curvature, const TargetT& current_image_estimate);
#if 0
  void compute_Hessian(TargetT& prior_Hessian_for_single_densel, 
                       const BasicCoordinate<3,int>& coords, const unsigned int input_param_num, 
                       const TargetT &current_image_estimate);
#endif

  virtual Succeeded add_multiplication_with_approximate_Hessian(TargetT& output, const TargetT& input) const;

  Array<3, float> get_weights() const;

  void set_weights(const Array<3, float>&);


  shared_ptr<TargetT> get_kappa_sptr() const;

  void set_kappa_sptr(const shared_ptr<TargetT>&);

  virtual Succeeded set_up(shared_ptr<const TargetT> const& target_sptr);

protected:
  bool only_2D;

  std::string gradient_filename_prefix;


  mutable Array<3, float> weights;
  std::string kappa_filename; // CHECK IF THERE IS A CONFILCT WHEN GIVING A KAPPA FILE...// THINK IF IT IS BETTER TO ESTIMATE
                              // KAPPA FILE IN THE CODE...

  virtual void check(TargetT const& current_image_estimate) const;

  virtual void set_defaults();
  virtual void initialise_keymap();
  virtual bool post_processing();
};

END_NAMESPACE_STIR

#endif