GibbsRelativeDifferencePenalty.h

Go to the documentation of this file.

//
//
/*
    Copyright (C) 2019 -2025, University College London
    Copyright (C) 2025, University of Milano-Bicocca
    This file is part of STIR.

    SPDX-License-Identifier: Apache-2.0

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

#include "stir/recon_buildblock/GibbsPenalty.h"
#include "stir/RegisteredParsingObject.h"
#include <cmath>
#include "stir/cuda_utilities.h"

#ifdef STIR_WITH_CUDA
#  include "stir/recon_buildblock/CUDA/CudaGibbsPenalty.h"
#endif

START_NAMESPACE_STIR

template <typename elemT>
class RelativeDifferencePotential
{
public:
  float gamma;
  float epsilon;
  __host__ __device__ inline double value(const elemT val_center, const elemT val_neigh, int z, int y, int x) const
  {
    // Implemented formula:
    // return 0.5 * (val_center -val_neigh)**2 / (val_center +val_neigh + gamma * |val_center - val_neigh| + epsilon);
    const elemT diff = val_center - val_neigh;
    const elemT add = val_center + val_neigh;
    const elemT NUM = 0.5 * (diff * diff);
    const elemT DEN = 1.0 / (add + gamma * fabs(diff) + epsilon);

    return NUM * DEN;
  }

  __host__ __device__ inline double derivative_10(const elemT val_center, const elemT val_neigh, int z, int y, int x) const
  {
    // Implemented formula:
    // return 0.5 * (val_center-val_neigh) * (val_center + 3*val_neigh + gamma * |val_center - val_neigh| + 2 * epsilon )  /
    // ((val_center +val_neigh)+ gamma * |val_center - val_neigh| + epsilon)**2;
    const elemT diff = val_center - val_neigh;
    const elemT factor = val_center + val_neigh + gamma * fabs(diff) + epsilon;
    const elemT NUM = 0.5 * diff * (factor + 2 * val_neigh + epsilon);
    const elemT DEN = 1.0 / (factor * factor);

    return NUM * DEN;
  }

  __host__ __device__ inline double derivative_20(const elemT val_center, const elemT val_neigh, int z, int y, int x) const
  {
    // Implemented formula:
    // return   (2*val_center + epsilon)**2 / (val_center + val_neigh + gamma * |val_center - val_neigh| + epsilon)**3;
    const elemT NUM = 2 * val_neigh + epsilon;
    const elemT DEN = 1.0 / (val_center + val_neigh + gamma * fabs(val_center - val_neigh) + epsilon);

    return NUM * NUM * DEN * DEN * DEN;
  }
  __host__ __device__ inline double derivative_11(const elemT val_center, const elemT val_neigh, int z, int y, int x) const
  {
    // Implemented formula:
    // return   -(2*val_center + epsilon)*(2*val_neigh + epsilon)/ (val_center + val_neigh + gamma * |val_center - val_neigh| +
    // epsilon)**3;
    const elemT NUM = -(2 * val_center + epsilon) * (2 * val_neigh + epsilon);
    const elemT DEN = 1.0 / (val_center + val_neigh + gamma * fabs(val_center - val_neigh) + epsilon);

    return NUM * DEN * DEN * DEN;
  }

  static inline bool is_convex() { return true; }

  void initialise_keymap(KeyParser& parser)
  {
    parser.add_key("gamma value", &this->gamma);
    parser.add_key("epsilon value", &this->epsilon);
  }

  void set_defaults()
  {
    this->gamma = 2;
    this->epsilon = 1e-7;
  }
};

template <typename elemT>
class GibbsRelativeDifferencePenalty : public RegisteredParsingObject<GibbsRelativeDifferencePenalty<elemT>,
                                                                      GeneralisedPrior<DiscretisedDensity<3, elemT>>,
                                                                      GibbsPenalty<elemT, RelativeDifferencePotential<elemT>>>
{
private:
  typedef RegisteredParsingObject<GibbsRelativeDifferencePenalty<elemT>,
                                  GeneralisedPrior<DiscretisedDensity<3, elemT>>,
                                  GibbsPenalty<elemT, RelativeDifferencePotential<elemT>>>
      base_type;

public:
  GibbsRelativeDifferencePenalty() { this->set_defaults(); }
  GibbsRelativeDifferencePenalty(const bool only_2D, float penalisation_factor, float gamma_v, float epsilon_v)
      : base_type(only_2D, penalisation_factor)
  {
    this->potential.gamma = gamma_v;
    this->potential.epsilon = epsilon_v;
  }

  static constexpr const char* registered_name = "Gibbs Relative Difference";
  float get_gamma() const { return this->potential.gamma; }
  float get_epsilon() const { return this->potential.epsilon; }
  void set_gamma(float gamma_v) { this->potential.gamma = gamma_v; }
  void set_epsilon(float epsilon_v) { this->potential.epsilon = epsilon_v; }
};

#ifdef STIR_WITH_CUDA

template <typename elemT>
class CudaGibbsRelativeDifferencePenalty
    : public RegisteredParsingObject<CudaGibbsRelativeDifferencePenalty<elemT>,
                                     GeneralisedPrior<DiscretisedDensity<3, elemT>>,
                                     CudaGibbsPenalty<elemT, RelativeDifferencePotential<elemT>>>
{
private:
  typedef RegisteredParsingObject<CudaGibbsRelativeDifferencePenalty<elemT>,
                                  GeneralisedPrior<DiscretisedDensity<3, elemT>>,
                                  CudaGibbsPenalty<elemT, RelativeDifferencePotential<elemT>>>
      base_type;

public:
  CudaGibbsRelativeDifferencePenalty() { this->set_defaults(); }
  CudaGibbsRelativeDifferencePenalty(const bool only_2D, float penalisation_factor, float gamma_v, float epsilon_v)
      : base_type(only_2D, penalisation_factor)
  {
    this->potential.gamma = gamma_v;
    this->potential.epsilon = epsilon_v;
  }

  static constexpr const char* registered_name = "Cuda Gibbs Relative Difference";
  float get_gamma() const { return this->potential.gamma; }
  float get_epsilon() const { return this->potential.epsilon; }
  void set_gamma(float gamma_v) { this->potential.gamma = gamma_v; }
  void set_epsilon(float epsilon_v) { this->potential.epsilon = epsilon_v; }
};
#endif

END_NAMESPACE_STIR

#endif