GibbsQuadraticPenalty.h

Go to the documentation of this file.

//
//
/*
    Copyright (C) 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_GibbsQuadraticPenalty_H__
#define __stir_recon_buildblock_GibbsQuadraticPenalty_H__

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

START_NAMESPACE_STIR

template <typename elemT>
class QuadraticPotential
{
public:
  __host__ __device__ inline double value(const elemT& val_center, const elemT val_neigh, int z, int y, int x) const
  {
    const elemT diff = val_center - val_neigh;
    return static_cast<double>(diff * diff) / 4.0;
  }
  __host__ __device__ inline double derivative_10(const elemT val_center, const elemT val_neigh, int z, int y, int x) const
  {
    return static_cast<double>(val_center - val_neigh) / 2.0;
  }
  __host__ __device__ inline double derivative_20(const elemT val_center, const elemT val_neigh, int z, int y, int x) const
  {
    return static_cast<double>(0.5);
  }
  __host__ __device__ inline double derivative_11(const elemT val_center, const elemT val_neigh, int z, int y, int x) const
  {
    return static_cast<double>(-0.5);
  }

  static inline bool is_convex() { return true; }
  void initialise_keymap(KeyParser& parser)
  {
    // No additional parameters needed for quadratic potential
  }
  void set_defaults()
  {
    // No additional parameters for quadratic potential
  }
};

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

public:
  static constexpr const char* const registered_name = "Gibbs Quadratic";

  GibbsQuadraticPenalty() { this->set_defaults(); }
  GibbsQuadraticPenalty(const bool only_2D, float penalisation_factor)
      : base_type(only_2D, penalisation_factor)
  {}
};

#ifdef STIR_WITH_CUDA

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

public:
  static constexpr const char* const registered_name = "Cuda Gibbs Quadratic";

  CudaGibbsQuadraticPenalty() { this->set_defaults(); }
  CudaGibbsQuadraticPenalty(const bool only_2D, float penalisation_factor)
      : base_type(only_2D, penalisation_factor)
  {}
};
#endif

END_NAMESPACE_STIR

#endif