LogcoshPrior.h

Go to the documentation of this file.

//
//
/*
 Copyright (C) 2000- 2011, Hammersmith Imanet Ltd
 Copyright (C) 2016, 2020, UCL
 This file is part of STIR.

 SPDX-License-Identifier: Apache-2.0

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

#include "stir/RegisteredParsingObject.h"
#include "stir/recon_buildblock/PriorWithParabolicSurrogate.h"
#include "stir/Array.h"
#include "stir/DiscretisedDensity.h"
#include "stir/shared_ptr.h"
#include <string>

START_NAMESPACE_STIR

template <typename elemT>
class LogcoshPrior : public RegisteredParsingObject<LogcoshPrior<elemT>,
                                                    GeneralisedPrior<DiscretisedDensity<3, elemT>>,
                                                    PriorWithParabolicSurrogate<DiscretisedDensity<3, elemT>>>
{
private:
  typedef RegisteredParsingObject<LogcoshPrior<elemT>,
                                  GeneralisedPrior<DiscretisedDensity<3, elemT>>,
                                  PriorWithParabolicSurrogate<DiscretisedDensity<3, elemT>>>
      base_type;

public:
  static const char* const registered_name;

  LogcoshPrior();

  LogcoshPrior(const bool only_2D, float penalization_factor);

  LogcoshPrior(const bool only_2D, float penalization_factor, const float scalar);

  bool parabolic_surrogate_curvature_depends_on_argument() const override { return false; }

  bool is_convex() const override;

  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;

  void parabolic_surrogate_curvature(DiscretisedDensity<3, elemT>& parabolic_surrogate_curvature,
                                     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 accumulate_Hessian_times_input(DiscretisedDensity<3, elemT>& output,
                                      const DiscretisedDensity<3, elemT>& current_estimate,
                                      const DiscretisedDensity<3, elemT>& input) const override;

  Array<3, float> get_weights() const;

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


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

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

  float get_scalar() const;

  void set_scalar(float scalar_v);

protected:
  bool only_2D;

  float scalar;


  std::string gradient_filename_prefix;


  mutable Array<3, float> weights;

  std::string kappa_filename;

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

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

private:
  shared_ptr<DiscretisedDensity<3, elemT>> kappa_ptr;


  static inline float logcosh(const float d)
  {
    const float x = fabs(d);
    if (x < 30.f)
      {
        return log(cosh(x));
      }
    else
      {
        return x + log(0.5f);
      }
  }


  static inline float surrogate(const float d, const float scalar)
  {
    const float eps = 0.01;
    const float x = d * scalar;
    // If abs(x) is less than eps,
    // use Taylor approximatation of tanh: tanh(x)/x ~= (x - x^3/3)/x = 1- x^2/3.
    // Prevents divide by zeros
    if (fabs(x) < eps)
      {
        return 1 - square(x) / 3;
      }
    else
      {
        return tanh(x) / x;
      }
  }


  elemT derivative_20(const elemT x_j, const elemT x_k) const;
  elemT derivative_11(const elemT x_j, const elemT x_k) const;
};

END_NAMESPACE_STIR

#endif