GibbsPenalty.h

Go to the documentation of this file.

/*
    Copyright (C) 2025, University College London
    Copyright (C) 2025, University of Milano-Bicocca
    Copyright (C) 2000- 2011, 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_GibbsPenalty_H__
#define __stir_recon_buildblock_GibbsPenalty_H__
 
#include "stir/recon_buildblock/GeneralisedPrior.h"
#include "stir/CartesianCoordinate3D.h"
#include "stir/Array.h"
#include "stir/DiscretisedDensity.h"
#include "stir/shared_ptr.h"
#include <string>
 
START_NAMESPACE_STIR

template <typename elemT, typename potentialT>
class GibbsPenalty : public GeneralisedPrior<DiscretisedDensity<3, elemT>>
{
private:
  typedef GeneralisedPrior<DiscretisedDensity<3, elemT>> base_type;
 
public:
  GibbsPenalty();

  GibbsPenalty(const bool only_2D, float penalization_factor);

  double compute_value(const DiscretisedDensity<3, elemT>& current_image_estimate) override;

  void compute_gradient(DiscretisedDensity<3, elemT>& prior_gradient,
                        const DiscretisedDensity<3, elemT>& current_image_estimate) override;

  double compute_gradient_times_input(const DiscretisedDensity<3, elemT>& input,
                                      const DiscretisedDensity<3, elemT>& current_image_estimate) override;

  void compute_Hessian(DiscretisedDensity<3, elemT>& prior_Hessian_for_single_densel,
                       const BasicCoordinate<3, int>& coords,
                       const DiscretisedDensity<3, elemT>& current_image_estimate) const override;

  void compute_Hessian_diagonal(DiscretisedDensity<3, elemT>& Hessian_diagonal,
                                const DiscretisedDensity<3, elemT>& current_estimate) const override;

  void accumulate_Hessian_times_input(DiscretisedDensity<3, elemT>& output,
                                      const DiscretisedDensity<3, elemT>& current_estimate,
                                      const DiscretisedDensity<3, elemT>& input) const override;

  const Array<3, float>& get_weights() const;

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


  shared_ptr<const DiscretisedDensity<3, elemT>> get_kappa_sptr() const;

  virtual void set_kappa_sptr(const shared_ptr<const DiscretisedDensity<3, elemT>>&);

  Succeeded set_up(shared_ptr<const DiscretisedDensity<3, elemT>> const& target_sptr) override;


  virtual std::string get_parsing_name() const;

  bool is_convex() const override;
 
protected:
  CartesianCoordinate3D<int> image_dim;         
  CartesianCoordinate3D<int> image_max_indices; 
  CartesianCoordinate3D<int> image_min_indices; 

  Array<3, float> weights;

  bool only_2D;


  std::string gradient_filename_prefix;

  std::string kappa_filename;

  potentialT potential;

  void compute_default_weights(const CartesianCoordinate3D<float>& grid_spacing, bool only_2D);

  void check(DiscretisedDensity<3, elemT> const& current_image_estimate) const override;
 
  void set_defaults() override;
  void initialise_keymap() override;
  bool post_processing() override;

  shared_ptr<const DiscretisedDensity<3, elemT>> kappa_ptr;
};
 
END_NAMESPACE_STIR
 
#include "stir/recon_buildblock/GibbsPenalty.inl"
 
#endif