diff --git a/src/buildblock/recon_array_functions.cxx b/src/buildblock/recon_array_functions.cxx index 5aab00180a..8d38e482ec 100644 --- a/src/buildblock/recon_array_functions.cxx +++ b/src/buildblock/recon_array_functions.cxx @@ -374,12 +374,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; } diff --git a/src/include/stir/recon_buildblock/LogcoshPrior.h b/src/include/stir/recon_buildblock/LogcoshPrior.h index 1593681ee9..d30de9eff6 100644 --- a/src/include/stir/recon_buildblock/LogcoshPrior.h +++ b/src/include/stir/recon_buildblock/LogcoshPrior.h @@ -125,6 +125,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 57dea783d6..2b0691266f 100644 --- a/src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h +++ b/src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h @@ -47,7 +47,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 4028169276..69321b2629 100644 --- a/src/include/stir/recon_buildblock/QuadraticPrior.h +++ b/src/include/stir/recon_buildblock/QuadraticPrior.h @@ -121,7 +121,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 f14742d7b1..055753d7d0 100644 --- a/src/recon_buildblock/LogcoshPrior.cxx +++ b/src/recon_buildblock/LogcoshPrior.cxx @@ -515,6 +515,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; @@ -577,6 +578,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 904b4b0d1b..7cf2baeae0 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -659,6 +659,78 @@ 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) + { + // 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 diff --git a/src/swig/stir.i b/src/swig/stir.i index 5a12471444..0253eea510 100644 --- a/src/swig/stir.i +++ b/src/swig/stir.i @@ -1698,6 +1698,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;