From 40f6b61d2a7982570f62e4aef43c78098da2072a Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Mon, 28 Sep 2020 15:38:31 +0100 Subject: [PATCH 1/7] Modify LL computation to be KL --- src/buildblock/recon_array_functions.cxx | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/buildblock/recon_array_functions.cxx b/src/buildblock/recon_array_functions.cxx index a966fe795a..4d0d8b2c44 100644 --- a/src/buildblock/recon_array_functions.cxx +++ b/src/buildblock/recon_array_functions.cxx @@ -382,12 +382,15 @@ void accumulate_loglikelihood(Viewgram& projection_data, // std::cerr << "Zero at " << r << ", " << b <<'\n'; const float new_estimate = max(estimated_projections[r][b], - projection_data[r][b]/max_quotient); - if (projection_data[r][b]<=small_value) - sub_result += - double(new_estimate); - else - sub_result += projection_data[r][b]*log(double(new_estimate)) - double(new_estimate); - } + projection_data[r][b]/max_quotient); + + if (projection_data[r][b] > 0.0) + { + sub_result += - (projection_data[r][b] * log(projection_data[r][b] / double(new_estimate)) + double(new_estimate) - projection_data[r][b]); + } else { + sub_result += - (double(new_estimate)); + } + } result += sub_result; } From 3adec826790fd0ff6ea80649094b4ccef21af070 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Thu, 26 Nov 2020 11:06:35 +0000 Subject: [PATCH 2/7] Logcosh fixes Prior to this OSMAPOSL was returning the error: FactoryRegistry: key LogCosh not found in current registry 5 possible values are: FilterRootPrior None PLS Quadratic Relative Difference Prior This commit fixes this and renames Logcosh to Logcosh Prior --- src/recon_buildblock/LogcoshPrior.cxx | 2 +- src/recon_buildblock/recon_buildblock_registries.cxx | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/recon_buildblock/LogcoshPrior.cxx b/src/recon_buildblock/LogcoshPrior.cxx index 99bf3a5480..3256978e00 100644 --- a/src/recon_buildblock/LogcoshPrior.cxx +++ b/src/recon_buildblock/LogcoshPrior.cxx @@ -125,7 +125,7 @@ LogcoshPrior::set_defaults() template <> const char * const LogcoshPrior::registered_name = - "Logcosh"; + "Logcosh Prior"; template LogcoshPrior::LogcoshPrior() diff --git a/src/recon_buildblock/recon_buildblock_registries.cxx b/src/recon_buildblock/recon_buildblock_registries.cxx index 9309ce3a95..ed1f212bdc 100644 --- a/src/recon_buildblock/recon_buildblock_registries.cxx +++ b/src/recon_buildblock/recon_buildblock_registries.cxx @@ -96,6 +96,7 @@ static FilterRootPrior >::RegisterIt dummy4; static QuadraticPrior::RegisterIt dummy5; static PLSPrior::RegisterIt dummyPLS; static RelativeDifferencePrior::RegisterIt dummyRelativeDifference; +static LogcoshPrior::RegisterIt dummyLogcosh; static ProjMatrixByBinUsingRayTracing::RegisterIt dummy11; static ProjMatrixByBinUsingInterpolation::RegisterIt dummy12; From 266a62e2760d2a947c9e2aca2ca30647ed3c1930 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Thu, 26 Nov 2020 11:21:05 +0000 Subject: [PATCH 3/7] Add Logcosh to users guide --- documentation/STIR-UsersGuide.tex | 42 ++++++++++++++++++++++++++++++- 1 file changed, 41 insertions(+), 1 deletion(-) diff --git a/documentation/STIR-UsersGuide.tex b/documentation/STIR-UsersGuide.tex index 8b6deb2fa1..1796543f74 100644 --- a/documentation/STIR-UsersGuide.tex +++ b/documentation/STIR-UsersGuide.tex @@ -4200,7 +4200,7 @@ \subsubsection{ \end{verbatim} -{ \subsubsubsection{Relative Difference} +{ \subsubsubsection{Relative Difference} \label{sec:priors:Relative_Difference} } This implements a Relative Difference Prior RDP, proposed by J. Nuyts, et.al., 2002. The gradient of the prior is computed as follows: \[g_r = \sum_{dr} w_{dr} * @@ -4244,6 +4244,45 @@ \subsubsection{ \end{verbatim} +{ \subsubsubsection{Logcosh} \label{sec:priors:Logcosh} +} +This implements a Logcosh. The gradient of the prior is computed as follows: + \[g_r = \sum_{dr} w_{dr} * \frac{1}{s} * \tanh(s*(\lambda_{r} - \lambda_{dr})) * + \kappa_r * \kappa_{r+dr}\] + + \noindent where $\lambda$ is the image where the gradient is computed + and $r$ and $dr$ are indices and the sum + is over the neighbourhood where the weights $w_{dr}$ are non-zero. + $s$ (a.k.a. scalar) controls the transition between the quadratic (smooth) and linear (edge-preserving) nature of the prior. + + The $\kappa$ image is the same as that of \ref{sec:priors:Quadratic}. + + By default, a 3x3 or 3x3x3 neigbourhood is used where the weights are set to + x-voxel\_size divided by the Euclidean distance between the points. + +{ \subsubsubsubsection{Parameters} +} + These are the keywords that can be used in addition to the ones listed in \ref{sec:priors}.. + \begin{verbatim} + prior type := Logcosh Prior + Logcosh Prior Parameters:= + ; next defaults to 0, set to 1 for 2D inverse Euclidean weights, 0 for 3D + only 2D:= 0 + ; next can be used to set weights explicitly. Needs to be a 3D array (of floats). + ' value of only_2D is ignored + ; following example uses 2D 'nearest neighbour' penalty + ; weights:={{{0,1,0},{1,0,1},{0,1,0}}} + ; use next parameter to specify an image with penalisation factors (a la Fessler) + ; see class documentation for more info + ; kappa filename:= + ; use next parameter to get gradient images at every subiteration + ; see class documentation + ; scalar is a parameter used to control edge preservation + ; scalar value := 1 + END Logcosh Prior Parameters:= + \end{verbatim} + + { \subsubsubsection{PLS} This implements the anatomical penalty function, Parallel Level Sets (PLS), proposed by Matthias J. Ehrhardt et. al [Ehr16]. Note that PLS becomes smoothed TV when an uniform anatomical image is provided. @@ -4280,6 +4319,7 @@ \subsubsection{ only_2D:=0 END PLS Prior Parameters:= \end{verbatim} + \subsubsection{ Selecting different projector pairs} \label{sec:projectorpairs} From 5237ef116c3837669576370979e1815f07f0559b Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Thu, 15 Oct 2020 11:21:27 +0100 Subject: [PATCH 4/7] Add method stir.ToPriorWithParabolicSurrogate() to python. Input is a GeneralisedPrior. This will give access to methods such as `parabolic_surrogate_curvature` and `add_multiplication_with_approximate_Hessian` from the objective function object. --- src/swig/stir.i | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/swig/stir.i b/src/swig/stir.i index db400915d7..2aa8cba073 100644 --- a/src/swig/stir.i +++ b/src/swig/stir.i @@ -1656,6 +1656,13 @@ stir::RegisteredParsingObject< stir::LogcoshPrior, stir::PriorWithParabolicSurrogate >; %template (LogcoshPrior3DFloat) stir::LogcoshPrior; +// Allows for access to parabolic surrogate type prior methods ( THIS IS A TEMPORARY FIX ) +%inline %{template stir::PriorWithParabolicSurrogate * ToPriorWithParabolicSurrogate(stir::GeneralisedPrior *b) { + return dynamic_cast*>(b); +} +%} +%template(ToPriorWithParabolicSurrogate) ToPriorWithParabolicSurrogate; + %template (Reconstruction3DFloat) stir::Reconstruction; //%template () stir::Reconstruction; %template (IterativeReconstruction3DFloat) stir::IterativeReconstruction; From 3b3cab119a64674cbf5b5436f2cb1996dca4bbf8 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Thu, 26 Nov 2020 17:44:38 +0000 Subject: [PATCH 5/7] Revert Logcosh registered_name to Logcosh and update docs [ci skip] --- documentation/STIR-UsersGuide.tex | 2 +- src/recon_buildblock/LogcoshPrior.cxx | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/documentation/STIR-UsersGuide.tex b/documentation/STIR-UsersGuide.tex index 1796543f74..6a5408d85a 100644 --- a/documentation/STIR-UsersGuide.tex +++ b/documentation/STIR-UsersGuide.tex @@ -4264,7 +4264,7 @@ \subsubsection{ } These are the keywords that can be used in addition to the ones listed in \ref{sec:priors}.. \begin{verbatim} - prior type := Logcosh Prior + prior type := Logcosh Logcosh Prior Parameters:= ; next defaults to 0, set to 1 for 2D inverse Euclidean weights, 0 for 3D only 2D:= 0 diff --git a/src/recon_buildblock/LogcoshPrior.cxx b/src/recon_buildblock/LogcoshPrior.cxx index 3256978e00..99bf3a5480 100644 --- a/src/recon_buildblock/LogcoshPrior.cxx +++ b/src/recon_buildblock/LogcoshPrior.cxx @@ -125,7 +125,7 @@ LogcoshPrior::set_defaults() template <> const char * const LogcoshPrior::registered_name = - "Logcosh Prior"; + "Logcosh"; template LogcoshPrior::LogcoshPrior() From e93d3a078ade14a865b4d54278a9d33b37af7181 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Thu, 10 Dec 2020 10:39:17 +0000 Subject: [PATCH 6/7] parabolic_surrogate_curvature_times_input methods for QP and Logcosh --- .../stir/recon_buildblock/LogcoshPrior.h | 5 ++ .../PriorWithParabolicSurrogate.h | 14 +++- .../stir/recon_buildblock/QuadraticPrior.h | 8 +- src/recon_buildblock/LogcoshPrior.cxx | 74 +++++++++++++++++++ src/recon_buildblock/QuadraticPrior.cxx | 73 ++++++++++++++++++ 5 files changed, 172 insertions(+), 2 deletions(-) diff --git a/src/include/stir/recon_buildblock/LogcoshPrior.h b/src/include/stir/recon_buildblock/LogcoshPrior.h index ab426cecdd..bfcc237e90 100644 --- a/src/include/stir/recon_buildblock/LogcoshPrior.h +++ b/src/include/stir/recon_buildblock/LogcoshPrior.h @@ -137,6 +137,11 @@ class LogcoshPrior: public const BasicCoordinate<3,int>& coords, const DiscretisedDensity<3,elemT> ¤t_image_estimate); + //! compute the parabolic surrogate curvature for the prior of the current image estimate and multiply by input image + void parabolic_surrogate_curvature_times_input(DiscretisedDensity<3,elemT>& output, + const DiscretisedDensity<3,elemT> ¤t_image_estimate, + const DiscretisedDensity<3,elemT> &input_image); + //! Compute the multiplication of the hessian of the prior multiplied by the input. virtual Succeeded accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& current_estimate, diff --git a/src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h b/src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h index 8786db6c3f..240163d52c 100644 --- a/src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h +++ b/src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h @@ -55,7 +55,19 @@ class PriorWithParabolicSurrogate: parabolic_surrogate_curvature(TargetT& parabolic_surrogate_curvature, const TargetT ¤t_estimate) = 0; - //! A function that allows skipping some computations if the curvature is independent of the \c current_estimate + //! This will calculate the parabolic surrogate curvature of the current image estimate multiplied by the input image + /*! + Function is comparable to that of accumulate_Hessian_times_input() but instead of the Hessian, we use the parabolic + surrogate curvature. + For each voxel (f_{j}) of the image, this method computes + \beta \sum_{i\in J} d_{i,j} f_{i} + where J is the neighbourhood about voxel j and d_{i,j} is the surrogate function between j and i. + */ + virtual void parabolic_surrogate_curvature_times_input(TargetT& output, + const TargetT& current_image_estimate, + const TargetT& input_image) = 0; + + //! A function that allows skipping some computations if the curvature is independent of the \c current_estimate /*! Defaults to return \c true, but can be overloaded by the derived class. */ virtual bool diff --git a/src/include/stir/recon_buildblock/QuadraticPrior.h b/src/include/stir/recon_buildblock/QuadraticPrior.h index 0f0262e227..23ddd829e0 100644 --- a/src/include/stir/recon_buildblock/QuadraticPrior.h +++ b/src/include/stir/recon_buildblock/QuadraticPrior.h @@ -129,7 +129,13 @@ class QuadraticPrior: public const BasicCoordinate<3,int>& coords, const DiscretisedDensity<3,elemT> ¤t_image_estimate); - //! Call accumulate_Hessian_times_input + //! compute the parabolic surrogate curvature for the prior of the current image estimate and multiply by input image + void parabolic_surrogate_curvature_times_input(DiscretisedDensity<3,elemT>& output, + const DiscretisedDensity<3,elemT> ¤t_image_estimate, + const DiscretisedDensity<3,elemT> &input_image); + + + //! Call accumulate_Hessian_times_input virtual Succeeded add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& input) const; diff --git a/src/recon_buildblock/LogcoshPrior.cxx b/src/recon_buildblock/LogcoshPrior.cxx index c05f4dd3ff..e56342407c 100644 --- a/src/recon_buildblock/LogcoshPrior.cxx +++ b/src/recon_buildblock/LogcoshPrior.cxx @@ -523,6 +523,7 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, // the only difference is that parabolic_surrogate_curvature uses input==1 assert( output.has_same_characteristics(input)); + output.fill(0.f); if (this->penalisation_factor==0) { return Succeeded::yes; @@ -585,6 +586,79 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, return Succeeded::yes; } +template +void +LogcoshPrior:: +parabolic_surrogate_curvature_times_input(DiscretisedDensity<3,elemT>& output, + const DiscretisedDensity<3,elemT> ¤t_image_estimate, + const DiscretisedDensity<3,elemT> &input_image) +{ + assert( output.has_same_characteristics(input)); + assert( output.has_same_characteristics(current_image_estimate)); + output.fill(0.f); + if (this->penalisation_factor==0) + { + return; + } + + DiscretisedDensityOnCartesianGrid<3,elemT>& output_cast = + dynamic_cast &>(output); + + if (weights.get_length() ==0) + { + compute_weights(weights, output_cast.get_grid_spacing(), this->only_2D); + } + + const bool do_kappa = !is_null_ptr(kappa_ptr); + + if (do_kappa && !kappa_ptr->has_same_characteristics(input_image)) + error("LogcoshPrior: kappa image has not the same index range as the reconstructed image\n"); + + const int min_z = output.get_min_index(); + const int max_z = output.get_max_index(); + for (int z=min_z; z<=max_z; z++) + { + const int min_dz = max(weights.get_min_index(), min_z-z); + const int max_dz = min(weights.get_max_index(), max_z-z); + + const int min_y = output[z].get_min_index(); + const int max_y = output[z].get_max_index(); + + for (int y=min_y;y<= max_y;y++) + { + const int min_dy = max(weights[0].get_min_index(), min_y-y); + const int max_dy = min(weights[0].get_max_index(), max_y-y); + + const int min_x = output[z][y].get_min_index(); + const int max_x = output[z][y].get_max_index(); + + for (int x=min_x;x<= max_x;x++) + { + const int min_dx = max(weights[0][0].get_min_index(), min_x-x); + const int max_dx = min(weights[0][0].get_max_index(), max_x-x); + + elemT result = 0; + for (int dz=min_dz;dz<=max_dz;++dz) + for (int dy=min_dy;dy<=max_dy;++dy) + for (int dx=min_dx;dx<=max_dx;++dx) + { + elemT voxel_diff= current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx]; + elemT current = weights[dz][dy][dx] * + surrogate(voxel_diff, this->scalar) * input_image[z+dz][y+dy][x+dx]; + + if (do_kappa) + current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; + + result += current; + } + + output[z][y][x] += result * this->penalisation_factor; + } + } + } + return; +} + # ifdef _MSC_VER // prevent warning message on reinstantiation, // note that we get a linking error if we don't have the explicit instantiation below diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index 68603f2962..1c9f8a9966 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -667,6 +667,79 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, return Succeeded::yes; } +template +void +QuadraticPrior:: +parabolic_surrogate_curvature_times_input(DiscretisedDensity<3,elemT>& output, + const DiscretisedDensity<3,elemT> ¤t_image_estimate, + const DiscretisedDensity<3,elemT> &input_image) +{ + assert( output.has_same_characteristics(input)); + assert( output.has_same_characteristics(current_image_estimate)); + output.fill(0.f); + if (this->penalisation_factor==0) + { + return; + } + + DiscretisedDensityOnCartesianGrid<3,elemT>& output_cast = + dynamic_cast &>(output); + + if (weights.get_length() ==0) + { + compute_weights(weights, output_cast.get_grid_spacing(), this->only_2D); + } + + const bool do_kappa = !is_null_ptr(kappa_ptr); + + if (do_kappa && !kappa_ptr->has_same_characteristics(input_image)) + error("LogcoshPrior: kappa image has not the same index range as the reconstructed image\n"); + + const int min_z = output.get_min_index(); + const int max_z = output.get_max_index(); + for (int z=min_z; z<=max_z; z++) + { + const int min_dz = max(weights.get_min_index(), min_z-z); + const int max_dz = min(weights.get_max_index(), max_z-z); + + const int min_y = output[z].get_min_index(); + const int max_y = output[z].get_max_index(); + + for (int y=min_y;y<= max_y;y++) + { + const int min_dy = max(weights[0].get_min_index(), min_y-y); + const int max_dy = min(weights[0].get_max_index(), max_y-y); + + const int min_x = output[z][y].get_min_index(); + const int max_x = output[z][y].get_max_index(); + + for (int x=min_x;x<= max_x;x++) + { + const int min_dx = max(weights[0][0].get_min_index(), min_x-x); + const int max_dx = min(weights[0][0].get_max_index(), max_x-x); + + elemT result = 0; + for (int dz=min_dz;dz<=max_dz;++dz) + for (int dy=min_dy;dy<=max_dy;++dy) + for (int dx=min_dx;dx<=max_dx;++dx) + { + elemT voxel_diff= current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx]; + // The parabolic surrogate curvature of the QP is 1 + elemT current = weights[dz][dy][dx] * input_image[z+dz][y+dy][x+dx]; + + if (do_kappa) + current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; + + result += current; + } + + output[z][y][x] += result * this->penalisation_factor; + } + } + } +} + + # ifdef _MSC_VER // prevent warning message on reinstantiation, // note that we get a linking error if we don't have the explicit instantiation below From 69694efbcdb7abce80e35c8678724379cd4e405f Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Thu, 10 Dec 2020 11:03:36 +0000 Subject: [PATCH 7/7] Remove redundant computation --- src/recon_buildblock/QuadraticPrior.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index 1c9f8a9966..2ece51be37 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -723,7 +723,6 @@ parabolic_surrogate_curvature_times_input(DiscretisedDensity<3,elemT>& output, for (int dy=min_dy;dy<=max_dy;++dy) for (int dx=min_dx;dx<=max_dx;++dx) { - elemT voxel_diff= current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx]; // The parabolic surrogate curvature of the QP is 1 elemT current = weights[dz][dy][dx] * input_image[z+dz][y+dy][x+dx];