Skip to content

Commit

Permalink
Remove epsilon from RDP. Use if statement on x_i == 0 and x_j == 0
Browse files Browse the repository at this point in the history
  • Loading branch information
@robbietuk
robbietuk committed Sep 25, 2020
1 parent ef435cc commit fc05869
Show file tree
Hide file tree
Showing 2 changed files with 19 additions and 34 deletions.
17 changes: 4 additions & 13 deletions src/include/stir/recon_buildblock/RelativeDifferencePrior.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,14 +52,14 @@ START_NAMESPACE_STIR
\f[
g_r = \sum_{dr} w_{dr} *
\frac{\left(\lambda_{j}-\lambda_{k}\right)\left(\gamma\left|\lambda_{j}-\lambda_{k}\right|+\lambda_{j}+3 \lambda_{k} + 2 \epsilon \right)}
{\left(\lambda_{j}+\lambda_{k}+\gamma\left|\lambda_{j}-\lambda_{k}\right| + \epsilon \right)^{2}} *
\frac{\left(\lambda_{j}-\lambda_{k}\right)\left(\gamma\left|\lambda_{j}-\lambda_{k}\right|+\lambda_{j}+3 \lambda_{k} \right)}
{\left(\lambda_{j}+\lambda_{k}+\gamma\left|\lambda_{j}-\lambda_{k}\right|\right)^{2}} *
\kappa_r * \kappa_{r+dr}
\f]
where \f$\lambda\f$ is the image and \f$r\f$ and \f$dr\f$ are indices and the sum
is over the neighbourhood where the weights \f$w_{dr}\f$ are non-zero. \f$\gamma\f$ is
a smoothing scalar term and the \f$\epsilon\f$ is a small positive value included to prevent division by zero.
a smoothing scalar term.
For more details, see: <em> J. Nuyts, D. Bequ&eacute;, P. Dupont, and L. Mortelmans,
“A Concave Prior Penalizing Relative Differences for Maximum-a-Posteriori Reconstruction in Emission Tomography,”
vol. 49, no. 1, pp. 56–60, 2002. </em>
Expand All @@ -83,7 +83,6 @@ START_NAMESPACE_STIR
; following example uses 2D 'nearest neighbour' penalty
; weights:={{{0,1,0},{1,0,1},{0,1,0}}}
; gamma value :=
; epsilon value :=
; see class documentation for more info
; use next parameter to specify an image with penalisation factors (a la Fessler)
; kappa filename:=
Expand Down Expand Up @@ -118,7 +117,7 @@ class RelativeDifferencePrior: public
RelativeDifferencePrior();

//! Constructs it explicitly
RelativeDifferencePrior(const bool only_2D, float penalization_factor, float gamma, float epsilon);
RelativeDifferencePrior(const bool only_2D, float penalization_factor, float gamma);

virtual bool
parabolic_surrogate_curvature_depends_on_argument() const
Expand All @@ -142,11 +141,6 @@ class RelativeDifferencePrior: public
//! set the gamma value used in the RDP
void set_gamma(float e);

//! get the epsilon value used in RDP
float get_epsilon() const;
//! set the epsilon value used in the RDP
void set_epsilon(float e);


//! get penalty weights for the neigbourhood
Array<3,float> get_weights() const;
Expand All @@ -171,9 +165,6 @@ class RelativeDifferencePrior: public
//! Create variable gamma for Relative Difference Penalty
float gamma;

//! Create variable epsilon for Relative Difference Penalty
float epsilon;

//! can be set during parsing to restrict the weights to the 2D case
bool only_2D;
//! filename prefix for outputing the gradient whenever compute_gradient() is called.
Expand Down
36 changes: 15 additions & 21 deletions src/recon_buildblock/RelativeDifferencePrior.cxx
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
//
//
//
/*
Copyright (C) 2000- 2019, Hammersmith Imanet Ltd
Expand Down Expand Up @@ -56,7 +56,6 @@ RelativeDifferencePrior<elemT>::initialise_keymap()
this->parser.add_key("weights", &weights);
this->parser.add_key("gradient filename prefix", &gradient_filename_prefix);
this->parser.add_key("gamma value", &this->gamma);
this->parser.add_key("epsilon value", &this->epsilon);
this->parser.add_stop_key("END Relative Difference Prior Parameters");
}

Expand Down Expand Up @@ -139,7 +138,6 @@ RelativeDifferencePrior<elemT>::set_defaults()
this->kappa_ptr.reset();
this->weights.recycle();
this->gamma = 1;
this->epsilon = 0.0001;
}

template <>
Expand Down Expand Up @@ -167,28 +165,14 @@ RelativeDifferencePrior<elemT>::
set_gamma(float g)
{ this->gamma = g; }

// Return the value of epsilon - a RDP parameter
template <typename elemT>
float
RelativeDifferencePrior<elemT>::
get_epsilon() const
{ return this->epsilon; }

// Set the value of epsilon - a RDP parameter
template <typename elemT>
void
RelativeDifferencePrior<elemT>::
set_epsilon(float g)
{ this->epsilon = g; }


template <typename elemT>
RelativeDifferencePrior<elemT>::RelativeDifferencePrior(const bool only_2D_v, float penalisation_factor_v, float gamma_v, float epsilon_v)
RelativeDifferencePrior<elemT>::RelativeDifferencePrior(const bool only_2D_v, float penalisation_factor_v, float gamma_v)
: only_2D(only_2D_v)
{
this->penalisation_factor = penalisation_factor_v;
this->gamma = gamma_v;
this->epsilon = epsilon_v;
}


Expand Down Expand Up @@ -313,10 +297,15 @@ compute_value(const DiscretisedDensity<3,elemT> &current_image_estimate)
for (int dx=min_dx;dx<=max_dx;++dx)
{
elemT current;
if (current_image_estimate[z][y][x] == 0.0 && current_image_estimate[z+dz][y+dy][x+dx] == 0.0){
current = 0.0;
} else {
current = weights[dz][dy][dx] * 0.5 *
(pow(current_image_estimate[z][y][x]-current_image_estimate[z+dz][y+dy][x+dx],2)/
(current_image_estimate[z][y][x]+current_image_estimate[z+dz][y+dy][x+dx]
+ this->gamma * abs(current_image_estimate[z][y][x]-current_image_estimate[z+dz][y+dy][x+dx]) + this->epsilon ));
+ this->gamma * abs(current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx])));
}

if (do_kappa)
current *=
(*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx];
Expand Down Expand Up @@ -388,12 +377,17 @@ compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient,
{

elemT current;
if (current_image_estimate[z][y][x] == 0.0 && current_image_estimate[z+dz][y+dy][x+dx] == 0.0){
current = 0.0;
} else {
current = weights[dz][dy][dx] *
(((current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx]) *
(this->gamma * abs(current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx]) +
current_image_estimate[z][y][x] + 3 * current_image_estimate[z+dz][y+dy][x+dx] + 2 * this->epsilon))/
current_image_estimate[z][y][x] + 3 * current_image_estimate[z+dz][y+dy][x+dx]))/
(square((current_image_estimate[z][y][x] + current_image_estimate[z+dz][y+dy][x+dx]) +
this->gamma * abs(current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx]) + this->epsilon)));
this->gamma * abs(current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx]))));
}

if (do_kappa)
current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx];

Expand Down

0 comments on commit fc05869

Please sign in to comment.