diff --git a/documentation/release_5.0.htm b/documentation/release_5.0.htm index 3818f7ca0a..edd390777b 100644 --- a/documentation/release_5.0.htm +++ b/documentation/release_5.0.htm @@ -143,6 +143,12 @@

Changed functionality

  • ROOT file I/O improvements. An entire entry's tree is no longer loaded by default and instead individual branches are loaded as needed. ROOT file event processing is now up to 4x faster. In addition, there is now an option to check energy window information (defaulting to on). Futhermore, added functionality to exclude true and unscattered event types from list mode processing.
    • +
  • + Poisson log-likelihood gradient methods now use actual_compute_subset_gradient_without_penalty + to properly handle the objective function gradient computation at the distributable level. + Removes the subtraction of the subset sensitivity, which lead to a inconsistency between log-likelihood function + and gradient when use_subset_sensitivities was false. +
    • diff --git a/recon_test_pack/OSSPS_test_PM_QP_subsens1.par b/recon_test_pack/OSSPS_test_PM_QP_subsens1.par new file mode 100644 index 0000000000..c9bde1a966 --- /dev/null +++ b/recon_test_pack/OSSPS_test_PM_QP_subsens1.par @@ -0,0 +1,73 @@ +OSSPSParameters := +; sample file for OSSPS +; parameters used here are for illustrative purposes only +; i.e. they are not recommended values + +objective function type:= PoissonLogLikelihoodWithLinearModelForMeanAndProjData +PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:= + +input file := Utahscat600k_ca_seg4.hs +zero end planes of segment 0:= 1 +; if disabled, defaults to maximum segment number in the file +maximum absolute segment number to process := 3 + +; change to STIR 2.x default for compatibility +use subset sensitivities:=1 +sensitivity filename:= RPTsens_seg3_PM.hv + +projector pair type := Matrix + Projector Pair Using Matrix Parameters := + Matrix type := Ray Tracing + Ray tracing matrix parameters := + number of rays in tangential direction to trace for each bin := 2 +; restrict to cylindrical fov := 0 + End Ray tracing matrix parameters := + End Projector Pair Using Matrix Parameters := + +; additive sinogram:=my_fake_randoms.hs +prior type := quadratic + Quadratic Prior Parameters:= + penalisation factor := 0.5 + ; next defaults to 0, set to 1 for 2D inverse Euclidean weights, 0 for 3D + only 2D:= 0 + END Quadratic Prior Parameters:= + +End PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:= + +output filename prefix := my_test_image_OSSPS_PM_QP_subsens +; iteration scheme + +number of subsets:= 4 +;start at subset:= 0 +;start at subiteration number := 1 +number of subiterations:= 8 +Save estimates at subiteration intervals:= 8 +;write update image := 0 +;report objective function values interval := 1 + +; if next is disabled, defaults to image full of 1s (but that's not good for OSSPS) +; in particular, make sure it has the correct scale +initial estimate:= test_image_PM_QP_6.hv +enforce initial positivity condition := 1 + +; here start OSSPS specific values + +; values to use for the 'precomputed denominator' +; specify either procomputed denomiator or normalisation type + ; use the following if you have it already (e.g. from previous run) + ; note: setting the value to 1 will use an images full of ones. + ; precomputed denominator := my_precomputed_denominator.hv + +; specify relaxation scheme +; lambda = relaxation_parameter/ (1+relaxation_gamma*(subiteration_num/num_subsets) +relaxation parameter := 2 +relaxation gamma:=.1 + + +END := + + + + + + diff --git a/recon_test_pack/run_tests.sh b/recon_test_pack/run_tests.sh index 9ea59ef442..f123bf52e0 100755 --- a/recon_test_pack/run_tests.sh +++ b/recon_test_pack/run_tests.sh @@ -245,6 +245,21 @@ echo There were problems here!; ThereWereErrors=1; fi +${MPIRUN} ${INSTALL_DIR}OSSPS OSSPS_test_PM_QP_subsens1.par 1> OSSPS_PM_QP.log 2> OSSPS_PM_QP_stderr.log + +echo '---- Comparing output of OSSPS subiter 8 using subset sensitivity (should be identical up to tolerance)' +echo Running ${INSTALL_DIR}compare_image +# relax test for the outer-rim voxels as these turn out to be more unstable than the internal ones +if ${INSTALL_DIR}compare_image -t 0.002 test_image_OSSPS_PM_QP_8.hv my_test_image_OSSPS_PM_QP_subsens_8.hv -a + ${INSTALL_DIR}compare_image -r 1 test_image_OSSPS_PM_QP_8.hv my_test_image_OSSPS_PM_QP_subsens_8.hv +then + echo ---- This test seems to be ok !; + else + echo There were problems here!; + ThereWereErrors=1; +fi + + echo echo ------------- tests on stir_math and correct_projdata --------- echo "first make up some randoms (just a projdata full of 1)" diff --git a/recon_test_pack/test_image_OSSPS_PM_QP_8.v b/recon_test_pack/test_image_OSSPS_PM_QP_8.v index 1380e0d12d..a0b76ce7c3 100644 Binary files a/recon_test_pack/test_image_OSSPS_PM_QP_8.v and b/recon_test_pack/test_image_OSSPS_PM_QP_8.v differ diff --git a/src/include/stir/recon_buildblock/GeneralisedPrior.h b/src/include/stir/recon_buildblock/GeneralisedPrior.h index caa4bf4dab..31a3249967 100644 --- a/src/include/stir/recon_buildblock/GeneralisedPrior.h +++ b/src/include/stir/recon_buildblock/GeneralisedPrior.h @@ -60,6 +60,17 @@ class GeneralisedPrior: virtual void compute_gradient(DataT& prior_gradient, const DataT ¤t_estimate) =0; + //! This should computes a single row of the Hessian + /*! Default implementation just call error(). This function needs to be overridden by the + derived class. + This method computes the Hessian of a particular voxel, indicated by \c coords, of the current image estimate. + The default method should overwrite the values in \c prior_Hessian_for_single_densel. + */ + virtual Succeeded + compute_Hessian(DataT& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DataT& current_image_estimate) const; + //! This should compute the multiplication of the Hessian with a vector and add it to \a output /*! Default implementation just call error(). This function needs to be overridden by the derived class. @@ -89,6 +100,9 @@ class GeneralisedPrior: virtual Succeeded set_up(shared_ptr const& target_sptr); + //! Returns the status of the _is_convex variable + bool get_is_convex() const; + protected: float penalisation_factor; //! sets value for penalisation factor @@ -102,6 +116,9 @@ class GeneralisedPrior: virtual void check(DataT const& current_estimate) const; bool _already_set_up; + + //! Variable to indicate that the prior is a convex function, should be set in defaults by the derived class + bool _is_convex; }; END_NAMESPACE_STIR diff --git a/src/include/stir/recon_buildblock/LogcoshPrior.h b/src/include/stir/recon_buildblock/LogcoshPrior.h index a3bed89478..127e63288f 100644 --- a/src/include/stir/recon_buildblock/LogcoshPrior.h +++ b/src/include/stir/recon_buildblock/LogcoshPrior.h @@ -125,9 +125,10 @@ class LogcoshPrior: public const DiscretisedDensity<3,elemT> ¤t_image_estimate); //! compute Hessian - void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, - const BasicCoordinate<3,int>& coords, - const DiscretisedDensity<3,elemT> ¤t_image_estimate); + virtual Succeeded + compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const; //! Compute the multiplication of the hessian of the prior multiplied by the input. virtual Succeeded accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, @@ -224,18 +225,19 @@ class LogcoshPrior: public { return tanh(x)/x; } } - //! The Hessian of log(cosh()) is sech^2(x) = (1/cosh(x))^2 + //! The second partial derivatives of the LogCosh Prior /*! - This function returns the hessian of the logcosh function - * @param d the difference between the ith and jth voxel. - * @param scalar is the logcosh scalar value controlling the priors transition between the quadratic and linear behaviour - * @return the second derivative of the log-cosh function + Diagonal refers to the second derivative w.r.t. x_j only (i.e. diagonal of the Hessian) + Off-diagonal refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal of the Hessian) + For LogCosh, the off diagonal is the negative of the diagonal. + * @param x_j is the target voxel. + * @param x_k is the voxel in the neighbourhood. + * @return the second order partial derivatives of the LogCosh Prior */ - static inline float Hessian(const float d, const float scalar) - { - const float x = d * scalar; - return square((1/ cosh(x))); - } + //@{ + float diagonal_second_derivative(const float x_j, const float x_k) const; + float off_diagonal_second_derivative(const float x_j, const float x_k) const; + //@} }; diff --git a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData.h b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData.h index b523a51f46..f93cd58190 100644 --- a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData.h +++ b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData.h @@ -70,10 +70,6 @@ public RegisteredParsingObject get_exam_info_uptr_for_target() const; @@ -170,6 +166,13 @@ public RegisteredParsingObject& arg) template void PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData:: -compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, - const TargetT ¤t_estimate, - const int subset_num) +actual_compute_subset_gradient_without_penalty(TargetT& gradient, + const TargetT ¤t_estimate, + const int subset_num, + const bool add_sensitivity) { if (subset_num<0 || subset_num>=this->get_num_subsets()) error("compute_sub_gradient_without_penalty subset_num out-of-range error"); @@ -452,9 +453,10 @@ compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, this->_single_frame_obj_funcs[frame_num]. - compute_sub_gradient_without_penalty_plus_sensitivity(dyn_gradient[frame_num], - dyn_image_estimate[frame_num], - subset_num); + actual_compute_subset_gradient_without_penalty(dyn_gradient[frame_num], + dyn_image_estimate[frame_num], + subset_num, + add_sensitivity); } this->_patlak_plot_sptr->multiply_dynamic_image_with_model_gradient(gradient, diff --git a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.h b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.h index 108127c6e8..5e8cc4596b 100644 --- a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.h +++ b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.h @@ -66,14 +66,13 @@ START_NAMESPACE_STIR the sensitivity because (if \f$r=0\f$) it is the total probability of detecting a count (in any bin) originating from \f$v\f$. - This class computes the gradient as a sum of these two terms. The - sensitivity has to be computed by the virtual function - \c add_subset_sensitivity(). The sum is computed by - \c compute_sub_gradient_without_penalty_plus_sensitivity(). - - The reason for this is that the sensitivity is data-independent, and - can be computed only once. See also - PoissonLogLikelihoodWithLinearModelForMeanAndListModeData. + This class computes the gradient directly, via \c compute_sub_gradient_without_penalty(). + This method is utilised by the \c OSSPS algorithm in STIR. + However, an additional method (\c compute_sub_gradient_without_penalty_plus_sensitivity()) + is provided that computes the sum of the subset gradient (without penalty) and the sensitivity. + This method is utilised by the \c OSMAPOSL algorithm. + + See also \c PoissonLogLikelihoodWithLinearModelForMeanAndListModeData. \par Relation with Kullback-Leibler distance @@ -120,38 +119,40 @@ public GeneralisedObjectiveFunction //PoissonLogLikelihoodWithLinearModelForMean(); - //! Implementation in terms of compute_sub_gradient_without_penalty_plus_sensitivity() - /*! \warning If separate subsensitivities are not used, we just subtract the total - sensitivity divided by the number of subsets. - This is fine for some algorithms as the sum over all the subsets is - equal to gradient of the objective function (without prior). - Other algorithms do not behave very stable under this approximation - however. So, currently setup() will return an error if - !subsets_are_approximately_balanced() and subset sensitivities - are not used. - - \see get_use_subset_sensitivities() - */ + //! Compute the subset gradient of the (unregularised) objective function + /*! + Implementation in terms of actual_compute_sub_gradient_without_penalty() + This function is used by OSSPS may be used by other gradient ascent/descent algorithms + + This computes + \f[ + {\partial L \over \partial \lambda_v} = + \sum_b P_{bv} ({y_b \over Y_b} - 1) + \f] + (see the class general documentation). + The sum will however be restricted to a subset. + */ virtual void compute_sub_gradient_without_penalty(TargetT& gradient, const TargetT ¤t_estimate, const int subset_num); - //! This should compute the gradient of the (unregularised) objective function plus the (sub)sensitivity - /*! - This function is used for instance by OSMAPOSL. - - This computes - \f[ {\partial L \over \partial \lambda_v} + P_v = - \sum_b P_{bv} {y_b \over Y_b} - \f] - (see the class general documentation). - The sum will however be restricted to a subset. - */ + //! This should compute the subset gradient of the (unregularised) objective function plus the subset sensitivity + /*! + Implementation in terms of actual_compute_sub_gradient_without_penalty(). + This function is used for instance by OSMAPOSL. + + This computes + \f[ {\partial L \over \partial \lambda_v} + P_v = + \sum_b P_{bv} {y_b \over Y_b} + \f] + (see the class general documentation). + The sum will however be restricted to a subset. + */ virtual void compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, const TargetT ¤t_estimate, - const int subset_num) =0; + const int subset_num); //! set-up sensitivity etc if possible /*! If \c recompute_sensitivity is \c false, we will try to @@ -264,6 +265,27 @@ public GeneralisedObjectiveFunction */ void compute_sensitivities(); + //! computes the subset gradient of the objective function without the penalty (optional: add subset sensitivity) + /*! + If \c add_sensitivity is \c true, this computes + \f[ {\partial L \over \partial \lambda_v} + P_v = + \sum_b P_{bv} {y_b \over Y_b} + \f] + (see the class general documentation). + The sum will however be restricted to a subset. + + However, if \c add_sensitivity is \c false, this function will instead compute only the gradient + \f[ + {\partial L \over \partial \lambda_v} = + \sum_b P_{bv} ({y_b \over Y_b} - 1) + \f] + */ + virtual void + actual_compute_subset_gradient_without_penalty(TargetT& gradient, + const TargetT ¤t_estimate, + const int subset_num, + const bool add_sensitivity) = 0; + //! Sets defaults for parsing /*! Resets \c sensitivity_filename, \c subset_sensitivity_filenames to empty, \c recompute_sensitivity to \c false, and \c use_subset_sensitivities to false. diff --git a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.h b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.h index 5b2664cb0c..27ae39bef6 100644 --- a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.h +++ b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.h @@ -79,10 +79,11 @@ public RegisteredParsingObject const& target_sptr) template void PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:: -compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, - const TargetT ¤t_estimate, - const int subset_num) +actual_compute_subset_gradient_without_penalty(TargetT& gradient, + const TargetT ¤t_estimate, + const int subset_num, + const bool add_sensitivity) { assert(subset_num>=0); assert(subset_numnum_subsets); @@ -435,9 +436,10 @@ compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, gated_gradient[gate_num].end_all(), 0.F); this->_single_gate_obj_funcs[gate_num]. - compute_sub_gradient_without_penalty_plus_sensitivity(gated_gradient[gate_num], - gated_image_estimate[gate_num], - subset_num); + actual_compute_subset_gradient_without_penalty(gated_gradient[gate_num], + gated_image_estimate[gate_num], + subset_num, + add_sensitivity); } // if(this->_motion_correction_type==-1) this->_reverse_motion_vectors.warp_image(gradient,gated_gradient) ; diff --git a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.h b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.h index c886bc2f94..1ab752b665 100644 --- a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.h +++ b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.h @@ -78,10 +78,6 @@ public PoissonLogLikelihoodWithLinearModelForMean virtual Succeeded set_up(shared_ptr const& target_sptr); - virtual - void compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, - const TargetT ¤t_estimate, - const int subset_num)=0; //! time frame definitions /*! \todo This is currently used to be able to compute the gradient for one time frame. However, it probably does not belong here. diff --git a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.h b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.h index 013b9bdca0..2eaeabe0ac 100644 --- a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.h +++ b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.h @@ -66,11 +66,17 @@ typedef RegisteredParsingObject(); - //! This should compute the gradient of the objective function at the \a current_estimate overwriting \a gradient - virtual - void compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, - const TargetT ¤t_estimate, - const int subset_num); + //! Computes the gradient of the objective function at the \a current_estimate overwriting \a gradient. + /*! + \warning If add_sensitivity = false and use_subset_sensitivities = false will return an error + because the gradient will not be correct. Try use_subset_sensitivities = true. + */ + virtual + void actual_compute_subset_gradient_without_penalty(TargetT& gradient, + const TargetT ¤t_estimate, + const int subset_num, + const bool add_sensitivity); + virtual TargetT * construct_target_ptr() const; int set_num_subsets(const int new_num_subsets); diff --git a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h index 23ff58c916..f13db8e5c6 100644 --- a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h +++ b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h @@ -212,10 +212,11 @@ public RegisteredParsingObject get_exam_info_uptr_for_target() const; @@ -272,11 +273,11 @@ public RegisteredParsingObject ¤t_image_estimate); //! compute Hessian - void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, - const BasicCoordinate<3,int>& coords, - const DiscretisedDensity<3,elemT> ¤t_image_estimate); + virtual Succeeded + compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const; //! Call accumulate_Hessian_times_input virtual Succeeded @@ -179,6 +180,20 @@ class QuadraticPrior: public virtual bool post_processing(); private: shared_ptr > kappa_ptr; + + //! The second partial derivatives of the Quadratic Prior + /*! + Diagonal refers to the second derivative w.r.t. x_j only (i.e. diagonal of the Hessian) + Off-diagonal refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal of the Hessian) + For the Quadratic Prior, the off diagonal is the negative of the diagonal. + * @param x_j is the target voxel. + * @param x_k is the voxel in the neighbourhood. + * @return the second order partial derivatives of the Quadratic Prior + */ + //@{ + float diagonal_second_derivative(const float x_j, const float x_k) const; + float off_diagonal_second_derivative(const float x_j, const float x_k) const; + //@} }; diff --git a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h index b4eaf4c74e..d8a0707c33 100644 --- a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h +++ b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h @@ -125,11 +125,20 @@ class RelativeDifferencePrior: public void compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, const DiscretisedDensity<3,elemT> ¤t_image_estimate); + //! compute Hessian + virtual Succeeded compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const; virtual Succeeded add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& input) const; + //! 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, + const DiscretisedDensity<3,elemT>& input) const; + //! get the gamma value used in RDP float get_gamma() const; //! set the gamma value used in the RDP @@ -193,6 +202,20 @@ class RelativeDifferencePrior: public virtual bool post_processing(); private: shared_ptr > kappa_ptr; + + //! The second partial derivatives of the Relative Difference Prior + /*! + Diagonal refers to the second derivative w.r.t. x_j only (i.e. diagonal of the Hessian) + Off-diagonal refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal of the Hessian) + See J. Nuyts, et al., 2002, Equation 7. + * @param x_j is the target voxel. + * @param x_k is the voxel in the neighbourhood. + * @return the second order partial derivatives of the Relative Difference Prior + */ + //@{ + float diagonal_second_derivative(const float x_j, const float x_k) const; + float off_diagonal_second_derivative(const float x_j, const float x_k) const; + //@} }; diff --git a/src/recon_buildblock/GeneralisedPrior.cxx b/src/recon_buildblock/GeneralisedPrior.cxx index 5751c0a6b5..5ad4d81222 100644 --- a/src/recon_buildblock/GeneralisedPrior.cxx +++ b/src/recon_buildblock/GeneralisedPrior.cxx @@ -40,6 +40,7 @@ GeneralisedPrior::set_defaults() { _already_set_up = false; this->penalisation_factor = 0; + this->_is_convex = false; } template @@ -51,6 +52,20 @@ set_up(shared_ptr const&) return Succeeded::yes; } +template +Succeeded +GeneralisedPrior:: +compute_Hessian(TargetT& output, + const BasicCoordinate<3,int>& coords, + const TargetT& current_image_estimate) const +{ + if (this->get_is_convex()) + error("GeneralisedPrior:\n compute_Hessian implementation is not overloaded by your convex prior."); + else + error("GeneralisedPrior:\n compute_Hessian is not implemented because the prior is not convex."); + return Succeeded::no; +} + template Succeeded GeneralisedPrior:: @@ -81,6 +96,11 @@ void GeneralisedPrior::check(TargetT const& current_estimate) const error("The prior should already be set-up, but it's not."); } +template +bool GeneralisedPrior::get_is_convex() const { + return this->_is_convex; +} + # ifdef _MSC_VER // prevent warning message on instantiation of abstract class # pragma warning(disable:4661) diff --git a/src/recon_buildblock/LogcoshPrior.cxx b/src/recon_buildblock/LogcoshPrior.cxx index f14742d7b1..78fe8036d3 100644 --- a/src/recon_buildblock/LogcoshPrior.cxx +++ b/src/recon_buildblock/LogcoshPrior.cxx @@ -108,6 +108,7 @@ void LogcoshPrior::set_defaults() { base_type::set_defaults(); + this->_is_convex = true; this->only_2D = false; this->scalar = 1.0; this->kappa_ptr.reset(); @@ -376,18 +377,21 @@ compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, } template -void +Succeeded LogcoshPrior:: compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, - const DiscretisedDensity<3,elemT> ¤t_image_estimate) + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const { assert( prior_Hessian_for_single_densel.has_same_characteristics(current_image_estimate)); prior_Hessian_for_single_densel.fill(0); if (this->penalisation_factor==0) { - return; + return Succeeded::yes; } + + this->check(current_image_estimate); + const DiscretisedDensityOnCartesianGrid<3,elemT>& current_image_cast = dynamic_cast< const DiscretisedDensityOnCartesianGrid<3,elemT> &>(current_image_estimate); @@ -421,20 +425,37 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, for (int dy=min_dy;dy<=max_dy;++dy) for (int dx=min_dx;dx<=max_dx;++dx) { - // sech^2(x * scalar); sech(x) = 1/cosh(x) - elemT voxel_diff= current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx]; - elemT current = weights[dz][dy][dx] * Hessian(voxel_diff, this->scalar); - - if (do_kappa) - current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; - - diagonal += current; - prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = -current*this->penalisation_factor; + elemT current = 0.0; + if (dz == 0 && dy == 0 && dx == 0) + { + // The j == k case (diagonal Hessian element), which is a sum over the neighbourhood. + for (int ddz=min_dz;ddz<=max_dz;++ddz) + for (int ddy=min_dy;ddy<=max_dy;++ddy) + for (int ddx=min_dx;ddx<=max_dx;++ddx) + { + elemT diagonal_current = weights[ddz][ddy][ddx] * + diagonal_second_derivative(current_image_estimate[z][y][x], + current_image_estimate[z+ddz][y+ddy][x+ddx]); + if (do_kappa) + diagonal_current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+ddz][y+ddy][x+ddx]; + current += diagonal_current; + } + } + else + { + // The j != k vases (off-diagonal Hessian elements) + current = weights[dz][dy][dx] * off_diagonal_second_derivative(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]; + } + prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = + current*this->penalisation_factor; } - prior_Hessian_for_single_densel[z][y][x]= diagonal * this->penalisation_factor; + return Succeeded::yes; } + template void LogcoshPrior::parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature, @@ -561,8 +582,18 @@ accumulate_Hessian_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_estimate[z][y][x] - current_estimate[z+dz][y+dy][x+dx]; - elemT current = weights[dz][dy][dx] * Hessian( voxel_diff, this->scalar) * input[z+dz][y+dy][x+dx]; + elemT current = weights[dz][dy][dx]; + if (dz == dy == dz == 0) { + // The j == k case + current *= diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; + } else { + current *= (diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + + off_diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * + input[z + dz][y + dy][x + dx]); + } if (do_kappa) current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; @@ -577,6 +608,22 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, return Succeeded::yes; } +template +float +LogcoshPrior:: +diagonal_second_derivative(const float x_j, const float x_k) const +{ + return square((1/ cosh((x_j - x_k) * scalar))); +} + +template +float +LogcoshPrior:: +off_diagonal_second_derivative(const float x_j, const float x_k) const +{ + return - diagonal_second_derivative(x_j, x_k); +} + # 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/PLSPrior.cxx b/src/recon_buildblock/PLSPrior.cxx index 0b782546fa..daa09ca30c 100644 --- a/src/recon_buildblock/PLSPrior.cxx +++ b/src/recon_buildblock/PLSPrior.cxx @@ -127,6 +127,7 @@ void PLSPrior::set_defaults() { base_type::set_defaults(); + this->_is_convex = true; this->only_2D = false; this->alpha=1; this->eta=1; diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.cxx index d63701a535..3851a50938 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.cxx @@ -358,29 +358,18 @@ compute_sub_gradient_without_penalty(TargetT& gradient, const TargetT ¤t_estimate, const int subset_num) { - if (subset_num<0 || subset_num>=this->get_num_subsets()) - error("compute_sub_gradient_without_penalty subset_num out-of-range error"); - - this-> - compute_sub_gradient_without_penalty_plus_sensitivity(gradient, - current_estimate, - subset_num); - // compute gradient -= sub_sensitivity - { - typename TargetT::full_iterator gradient_iter = - gradient.begin_all(); - const typename TargetT::full_iterator gradient_end = - gradient.end_all(); - typename TargetT::const_full_iterator sensitivity_iter = - this->get_subset_sensitivity(subset_num).begin_all_const(); - while (gradient_iter != gradient_end) - { - *gradient_iter -= (*sensitivity_iter); - ++gradient_iter; ++sensitivity_iter; - } - } + this->actual_compute_subset_gradient_without_penalty(gradient, current_estimate, subset_num, false); } +template +void +PoissonLogLikelihoodWithLinearModelForMean:: +compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, + const TargetT ¤t_estimate, + const int subset_num) +{ + actual_compute_subset_gradient_without_penalty(gradient, current_estimate, subset_num, true); +} template void diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx index f59e0a77c3..d585cfb8d0 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx @@ -421,16 +421,20 @@ construct_target_ptr() const this->target_parameter_parser.create(this->get_input_data()); } -template -void -PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin:: -compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, - const TargetT ¤t_estimate, - const int subset_num) -{ - +template +void +PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin:: +actual_compute_subset_gradient_without_penalty(TargetT& gradient, + const TargetT ¤t_estimate, + const int subset_num, + const bool add_sensitivity) +{ assert(subset_num>=0); assert(subset_numnum_subsets); + if (!add_sensitivity && !this->get_use_subset_sensitivities()) + error("PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin::" + "actual_compute_subset_gradient_without_penalty(): cannot subtract subset sensitivity because " + "use_subset_sensitivities is false. This will result in an error in the gradient computation."); const double start_time = this->frame_defs.get_start_time(this->current_frame_num); const double end_time = this->frame_defs.get_end_time(this->current_frame_num); @@ -532,6 +536,23 @@ compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, } } info(boost::format("Number of used events: %1%") % num_used_events); + + if (!add_sensitivity) + { + // subtract the subset sensitivity + // compute gradient -= sub_sensitivity + typename TargetT::full_iterator gradient_iter = + gradient.begin_all(); + const typename TargetT::full_iterator gradient_end = + gradient.end_all(); + typename TargetT::const_full_iterator sensitivity_iter = + this->get_subset_sensitivity(subset_num).begin_all_const(); + while (gradient_iter != gradient_end) + { + *gradient_iter -= (*sensitivity_iter); + ++gradient_iter; ++sensitivity_iter; + } + } } # ifdef _MSC_VER diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx index 84bef7219f..e0591a55e0 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx @@ -610,9 +610,10 @@ set_up_before_sensitivity(shared_ptr const& target_sptr) template void PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: -compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, - const TargetT ¤t_estimate, - const int subset_num) +actual_compute_subset_gradient_without_penalty(TargetT& gradient, + const TargetT ¤t_estimate, + const int subset_num, + const bool add_sensitivity) { assert(subset_num>=0); assert(subset_numnum_subsets); @@ -625,14 +626,13 @@ compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, subset_num, this->num_subsets, -this->max_segment_num_to_process, - this->max_segment_num_to_process, - this->zero_seg0_end_planes!=0, - NULL, - this->additive_proj_data_sptr - , caching_info_ptr - ); - - + this->max_segment_num_to_process, + this->zero_seg0_end_planes!=0, + NULL, + this->additive_proj_data_sptr, + this->normalisation_sptr, + caching_info_ptr, + add_sensitivity); } @@ -1052,10 +1052,10 @@ actual_accumulate_sub_Hessian_times_input_without_penalty(TargetT& output, shared_ptr tmp(output.get_empty_copy()); this->get_projector_pair().get_back_projector_sptr()->get_output(*tmp); - // output += tmp; + // output -= tmp; std::transform(output.begin_all(), output.end_all(), tmp->begin_all(), output.begin_all(), - std::plus()); + std::minus()); return Succeeded::yes; } @@ -1077,7 +1077,7 @@ RPC_process_related_viewgrams_type RPC_process_related_viewgrams_sensitivity_com #else //! Call-back function for compute_gradient -static RPC_process_related_viewgrams_type RPC_process_related_viewgrams_gradient; +template static RPC_process_related_viewgrams_type RPC_process_related_viewgrams_gradient; //! Call-back function for accumulate_loglikelihood static RPC_process_related_viewgrams_type RPC_process_related_viewgrams_accumulate_loglikelihood; @@ -1097,10 +1097,13 @@ void distributable_compute_gradient(const shared_ptr& for bool zero_seg0_end_planes, double* log_likelihood_ptr, shared_ptr const& additive_binwise_correction, - DistributedCachingInformation* caching_info_ptr + shared_ptr const& normalisation_sptr, + DistributedCachingInformation* caching_info_ptr, + const bool add_sensitivity ) { - + if (add_sensitivity){ + // Within the RPC process, subtract ones before to back projection ( backproj[ y/ybar - 1] ) distributable_computation(forward_projector_sptr, back_projector_sptr, symmetries_sptr, @@ -1111,10 +1114,28 @@ void distributable_compute_gradient(const shared_ptr& for zero_seg0_end_planes, log_likelihood_ptr, additive_binwise_correction, - /* normalisation info to be ignored */ shared_ptr(), 0., 0., - &RPC_process_related_viewgrams_gradient, + /* normalisation info to be ignored */ shared_ptr(), + 0., 0., + &RPC_process_related_viewgrams_gradient, caching_info_ptr - ); + ); + } else if (!add_sensitivity){ + // Within the RPC process, only do div/truncate ( backproj[ y/ybar ] ) + distributable_computation(forward_projector_sptr, + back_projector_sptr, + symmetries_sptr, + &output_image, &input_image, + proj_dat, true, //i.e. do read projection data + subset_num, num_subsets, + min_segment, max_segment, + zero_seg0_end_planes, + log_likelihood_ptr, + additive_binwise_correction, + normalisation_sptr, 0., 0., + &RPC_process_related_viewgrams_gradient, + caching_info_ptr + ); + } } @@ -1195,7 +1216,7 @@ void distributable_sensitivity_computation( //////////// RPC functions - +template void RPC_process_related_viewgrams_gradient( const shared_ptr& forward_projector_sptr, const shared_ptr& back_projector_sptr, @@ -1205,8 +1226,6 @@ void RPC_process_related_viewgrams_gradient( const RelatedViewgrams* mult_viewgrams_ptr) { assert(measured_viewgrams_ptr != NULL); - if (!is_null_ptr(mult_viewgrams_ptr)) - error("Internal error: mult_viewgrams_ptr should be zero when computing gradient"); RelatedViewgrams estimated_viewgrams = measured_viewgrams_ptr->get_empty_copy(); @@ -1229,23 +1248,28 @@ void RPC_process_related_viewgrams_gradient( } */ forward_projector_sptr->forward_project(estimated_viewgrams); - - - + if (additive_binwise_correction_ptr != NULL) - { estimated_viewgrams += (*additive_binwise_correction_ptr); - } - - - - - // for sinogram division - divide_and_truncate(*measured_viewgrams_ptr, estimated_viewgrams, rim_truncation_sino, count, count2, log_likelihood_ptr); - + + // adding the sensitivity: backproj[y/ybar] * + // not adding the sensitivity computes the gradient: backproj[y/ybar - 1] * + // * ignoring normalisation * + if (!add_sensitivity){ + if (mult_viewgrams_ptr) + { + // subtract normalised ones from the data [y/ybar - 1/N] + *measured_viewgrams_ptr -= *mult_viewgrams_ptr; + } else { + // No mult_viewgrams_ptr, subtract ones [y/ybar - 1] + *measured_viewgrams_ptr -= 1; + } + } + + // back project back_projector_sptr->back_project(*measured_viewgrams_ptr); }; diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index 904b4b0d1b..d206f02071 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -123,6 +123,7 @@ void QuadraticPrior::set_defaults() { base_type::set_defaults(); + this->_is_convex = true; this->only_2D = false; this->kappa_ptr.reset(); this->weights.recycle(); @@ -430,17 +431,17 @@ compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, } template -void +Succeeded QuadraticPrior:: compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, - const DiscretisedDensity<3,elemT> ¤t_image_estimate) + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const { assert( prior_Hessian_for_single_densel.has_same_characteristics(current_image_estimate)); prior_Hessian_for_single_densel.fill(0); if (this->penalisation_factor==0) { - return; + return Succeeded::yes; } this->check(current_image_estimate); @@ -479,19 +480,34 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, for (int dy=min_dy;dy<=max_dy;++dy) for (int dx=min_dx;dx<=max_dx;++dx) { - // dz==0,dy==0,dx==0 will have weight 0, so we can just include it in the loop - elemT current = - weights[dz][dy][dx]; - - if (do_kappa) - current *= - (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; - - diagonal += current; - prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = -current*this->penalisation_factor; + elemT current = 0.0; + if (dz == 0 && dy == 0 && dx == 0) + { + // The j == k case (diagonal Hessian element), which is a sum over the neighbourhood. + for (int ddz=min_dz;ddz<=max_dz;++ddz) + for (int ddy=min_dy;ddy<=max_dy;++ddy) + for (int ddx=min_dx;ddx<=max_dx;++ddx) + { + elemT diagonal_current = weights[ddz][ddy][ddx] * + diagonal_second_derivative(current_image_estimate[z][y][x], + current_image_estimate[z+ddz][y+ddy][x+ddx]); + if (do_kappa) + diagonal_current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+ddz][y+ddy][x+ddx]; + current += diagonal_current; + } + } + else + { + // The j != k vases (off-diagonal Hessian elements) + current = weights[dz][dy][dx] * off_diagonal_second_derivative(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]; + } + prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = + current*this->penalisation_factor; } - prior_Hessian_for_single_densel[z][y][x]= diagonal * this->penalisation_factor; + return Succeeded::yes; } template @@ -581,14 +597,79 @@ QuadraticPrior:: add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& input) const { - return accumulate_Hessian_times_input(output, input, input); + // TODO this function overlaps enormously with parabolic_surrogate_curvature + // the only difference is that parabolic_surrogate_curvature uses input==1 + + assert( output.has_same_characteristics(input)); + if (this->penalisation_factor==0) + { + return Succeeded::yes; + } + + this->check(input); + + 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)) + error("QuadraticPrior: 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 current = + weights[dz][dy][dx] * input[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 Succeeded::yes; } template Succeeded QuadraticPrior:: accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, - const DiscretisedDensity<3,elemT>& /*current_estimate*/, + const DiscretisedDensity<3,elemT>& current_estimate, const DiscretisedDensity<3,elemT>& input) const { // TODO this function overlaps enormously with parabolic_surrogate_curvature @@ -638,13 +719,31 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, 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); + /// At this point, we have j = [z][y][x] + // The next for loops will have k = [z+dz][y+dy][x+dx] + // The following computes + //(H_{wf} y)_j = + // \sum_{k\in N_j} w_{(j,k)} f''_{d}(x_j,x_k) y_j + + // \sum_{(i \in N_j) \ne j} w_{(j,i)} f''_{od}(x_j, x_i) y_i + // Note the condition in the second sum that i is not equal to j + 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 current = - weights[dz][dy][dx] * input[z+dz][y+dy][x+dx]; + elemT current = weights[dz][dy][dx]; + if (dz == dy == dz == 0) { + // The j == k case + current *= diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; + } else { + current *= (diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + + off_diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * + input[z + dz][y + dy][x + dx]); + } if (do_kappa) current *= @@ -659,6 +758,18 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, return Succeeded::yes; } +template +float +QuadraticPrior:: +diagonal_second_derivative(const float x_j, const float x_k) const +{ return 1.0;} + +template +float +QuadraticPrior:: +off_diagonal_second_derivative(const float x_j, const float x_k) const +{ return -1.0;} + # 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/RelativeDifferencePrior.cxx b/src/recon_buildblock/RelativeDifferencePrior.cxx index e7fce65ea5..1fd446156b 100644 --- a/src/recon_buildblock/RelativeDifferencePrior.cxx +++ b/src/recon_buildblock/RelativeDifferencePrior.cxx @@ -127,6 +127,7 @@ void RelativeDifferencePrior::set_defaults() { base_type::set_defaults(); + this->_is_convex = true; this->only_2D = false; this->kappa_ptr.reset(); this->weights.recycle(); @@ -178,6 +179,7 @@ template RelativeDifferencePrior::RelativeDifferencePrior(const bool only_2D_v, float penalisation_factor_v, float gamma_v, float epsilon_v) : only_2D(only_2D_v) { + set_defaults(); this->penalisation_factor = penalisation_factor_v; this->gamma = gamma_v; this->epsilon = epsilon_v; @@ -420,6 +422,86 @@ compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, } } +template +Succeeded +RelativeDifferencePrior:: +compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const +{ + assert( prior_Hessian_for_single_densel.has_same_characteristics(current_image_estimate)); + prior_Hessian_for_single_densel.fill(0); + if (this->penalisation_factor==0) + { + return Succeeded::yes; + } + + this->check(current_image_estimate); + + const DiscretisedDensityOnCartesianGrid<3,elemT>& current_image_cast = + dynamic_cast< const DiscretisedDensityOnCartesianGrid<3,elemT> &>(current_image_estimate); + + DiscretisedDensityOnCartesianGrid<3,elemT>& prior_Hessian_for_single_densel_cast = + dynamic_cast &>(prior_Hessian_for_single_densel); + + if (weights.get_length() ==0) + { + compute_weights(weights, current_image_cast.get_grid_spacing(), this->only_2D); + } + + + const bool do_kappa = !is_null_ptr(kappa_ptr); + + if (do_kappa && kappa_ptr->has_same_characteristics(current_image_estimate)) + error("RelativeDifferencePrior: kappa image has not the same index range as the reconstructed image\n"); + + const int z = coords[1]; + const int y = coords[2]; + const int x = coords[3]; + const int min_dz = max(weights.get_min_index(), prior_Hessian_for_single_densel.get_min_index()-z); + const int max_dz = min(weights.get_max_index(), prior_Hessian_for_single_densel.get_max_index()-z); + + const int min_dy = max(weights[0].get_min_index(), prior_Hessian_for_single_densel[z].get_min_index()-y); + const int max_dy = min(weights[0].get_max_index(), prior_Hessian_for_single_densel[z].get_max_index()-y); + + const int min_dx = max(weights[0][0].get_min_index(), prior_Hessian_for_single_densel[z][y].get_min_index()-x); + const int max_dx = min(weights[0][0].get_max_index(), prior_Hessian_for_single_densel[z][y].get_max_index()-x); + + elemT diagonal = 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 current = 0.0; + if (dz == 0 && dy == 0 && dx == 0) + { + // The j == k case (diagonal Hessian element), which is a sum over the neighbourhood. + for (int ddz=min_dz;ddz<=max_dz;++ddz) + for (int ddy=min_dy;ddy<=max_dy;++ddy) + for (int ddx=min_dx;ddx<=max_dx;++ddx) + { + elemT diagonal_current = weights[ddz][ddy][ddx] * + diagonal_second_derivative(current_image_estimate[z][y][x], + current_image_estimate[z+ddz][y+ddy][x+ddx]); + if (do_kappa) + diagonal_current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+ddz][y+ddy][x+ddx]; + current += diagonal_current; + } + } + else + { + // The j != k vases (off-diagonal Hessian elements) + current = weights[dz][dy][dx] * off_diagonal_second_derivative(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]; + } + prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = + current*this->penalisation_factor; + } + + return Succeeded::yes; +} + template Succeeded RelativeDifferencePrior:: @@ -430,6 +512,121 @@ add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, return Succeeded::no; } +template +Succeeded +RelativeDifferencePrior:: +accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, + const DiscretisedDensity<3,elemT>& current_estimate, + const DiscretisedDensity<3,elemT>& input) const +{ + // TODO this function overlaps enormously with parabolic_surrogate_curvature + // the only difference is that parabolic_surrogate_curvature uses input==1 + + assert( output.has_same_characteristics(input)); + if (this->penalisation_factor==0) + { + return Succeeded::yes; + } + + 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)) + 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); + + /// At this point, we have j = [z][y][x] + // The next for loops will have k = [z+dz][y+dy][x+dx] + // The following computes + //(H_{wf} y)_j = + // \sum_{k\in N_j} w_{(j,k)} f''_{d}(x_j,x_k) y_j + + // \sum_{(i \in N_j) \ne j} w_{(j,i)} f''_{od}(x_j, x_i) y_i + // Note the condition in the second sum that i is not equal to j + + 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 current = weights[dz][dy][dx]; + if (dz == dy == dz == 0) { + // The j == k case + current *= diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; + } else { + current *= (diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + + off_diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * + input[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 Succeeded::yes; +} + +template +float +RelativeDifferencePrior:: +diagonal_second_derivative(const float x_j, const float x_k) const +{ + if (x_j > 0.0 || x_k > 0.0 || this->epsilon > 0.0) + return 2 * pow(2 * x_k + this->epsilon, 2) / + pow(x_j + x_k + this->gamma * abs(x_j - x_k) + this->epsilon, 3); + else + return 0.0; +} + +template +float +RelativeDifferencePrior:: +off_diagonal_second_derivative(const float x_j, const float x_k) const +{ + if (x_j > 0.0 || x_k > 0.0 || this->epsilon > 0.0) + return - 2 * (2 * x_j + this->epsilon)*(2 * x_k + this->epsilon) / + pow(x_j + x_k + this->gamma * abs(x_j - x_k) + this->epsilon, 3); + else + return 0.0; +} + # 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/distributable.cxx b/src/recon_buildblock/distributable.cxx index a7f034737e..129ce51d99 100644 --- a/src/recon_buildblock/distributable.cxx +++ b/src/recon_buildblock/distributable.cxx @@ -216,6 +216,13 @@ void get_viewgrams(shared_ptr >& y, #endif normalisation_sptr->undo(*mult_viewgrams_sptr); } + else if (zero_seg0_end_planes) + { + // No normalisation provided but zero_seg0_end_planes, create a mult_viewgrams + mult_viewgrams_sptr.reset( + new RelatedViewgrams(proj_dat_ptr->get_empty_related_viewgrams(view_segment_num, symmetries_ptr))); + mult_viewgrams_sptr->fill(1.F); + } if (view_segment_num.segment_num()==0 && zero_seg0_end_planes) { diff --git a/src/recon_buildblock/distributableMPICacheEnabled.cxx b/src/recon_buildblock/distributableMPICacheEnabled.cxx index 0bea637b5e..cd65441907 100644 --- a/src/recon_buildblock/distributableMPICacheEnabled.cxx +++ b/src/recon_buildblock/distributableMPICacheEnabled.cxx @@ -126,6 +126,13 @@ void get_viewgrams(shared_ptr >& y, mult_viewgrams_sptr->fill(1.F); normalisation_sptr->undo(*mult_viewgrams_sptr); } + else if (zero_seg0_end_planes) + { + // No normalisation provided but zero_seg0_end_planes, create a mult_viewgrams + mult_viewgrams_sptr.reset( + new RelatedViewgrams(proj_dat_ptr->get_empty_related_viewgrams(view_segment_num, symmetries_ptr))); + mult_viewgrams_sptr->fill(1.F); + } if (view_segment_num.segment_num()==0 && zero_seg0_end_planes) { diff --git a/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx b/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx index 7d43a549f9..4dea600568 100644 --- a/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx +++ b/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx @@ -107,6 +107,14 @@ class PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests : public RunTes /*! Note that this function is not specific to PoissonLogLikelihoodWithLinearModelForMeanAndProjData */ void run_tests_for_objective_function(GeneralisedObjectiveFunction& objective_function, target_type& target); + + //! Test the gradient of the objective function by comparing to the numerical gradient via perturbation + void test_objective_function_gradient(GeneralisedObjectiveFunction& objective_function, + target_type& target); + + //! Test the Hessian of the objective function by testing the (x^T Hx > 0) condition + void test_objective_function_Hessian_concavity(GeneralisedObjectiveFunction& objective_function, + target_type& target); }; PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests:: @@ -118,6 +126,20 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests:: run_tests_for_objective_function(GeneralisedObjectiveFunction& objective_function, PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::target_type& target) { + std::cerr << "----- testing Gradient\n"; + test_objective_function_gradient(objective_function, + target); + + std::cerr << "----- testing Hessian-vector product (accumulate_Hessian_times_input)\n"; + test_objective_function_Hessian_concavity(objective_function, + target); +} + +void +PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests:: +test_objective_function_gradient(GeneralisedObjectiveFunction &objective_function, + target_type &target) +{ shared_ptr gradient_sptr(target.get_empty_copy()); shared_ptr gradient_2_sptr(target.get_empty_copy()); const int subset_num = 0; @@ -166,6 +188,38 @@ run_tests_for_objective_function(GeneralisedObjectiveFunction &objective_function, + target_type &target){ + /// setup images + shared_ptr output(target.get_empty_copy()); + + /// Compute H x + objective_function.accumulate_Hessian_times_input(*output, target, target); + + /// Compute dot(x,(H x)) + float my_sum = 0.0; + { + target_type::full_iterator target_iter = target.begin_all(); + target_type::full_iterator output_iter = output->begin_all(); + while (target_iter != target.end_all())// && testOK) + { + my_sum += *target_iter * *output_iter; + ++target_iter; ++output_iter; + } + } + // test for a CONCAVE function + if (this->check_if_less( my_sum, 0)) { +// info("PASS: Computation of x^T H x = " + std::to_string(my_sum) + " < 0" and is therefore concave); + } else { + // print to console the FAILED configuration + info("FAIL: Computation of x^T H x = " + std::to_string(my_sum) + " > 0 and is therefore NOT concave" + + "\n >target image max=" + std::to_string(target.find_max()) + + "\n >target image min=" + std::to_string(target.find_min())); + } +} + void PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests:: diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index 961afdf0a8..b3f7d0a1e7 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -29,11 +29,12 @@ #include "stir/recon_buildblock/QuadraticPrior.h" #include "stir/recon_buildblock/RelativeDifferencePrior.h" #include "stir/recon_buildblock/LogcoshPrior.h" -//#include "stir/recon_buildblock/PLSPrior.h" +#include "stir/recon_buildblock/PLSPrior.h" #include "stir/RunTests.h" #include "stir/IO/read_from_file.h" #include "stir/IO/write_to_file.h" #include "stir/info.h" +#include "stir/Verbosity.h" #include "stir/Succeeded.h" #include "stir/num_threads.h" #include @@ -53,6 +54,8 @@ START_NAMESPACE_STIR This test compares the result of GeneralisedPrior::compute_gradient() with a numerical gradient computed by using the GeneralisedPrior::compute_value() function. + Additionally, the Hessian's convexity is tested, via GeneralisedPrior::accumulate_Hessian_times_input(), + by evaluating the x^T Hx > 0 constraint. */ class GeneralisedPriorTests : public RunTests @@ -70,6 +73,14 @@ class GeneralisedPriorTests : public RunTests void construct_input_data(shared_ptr& density_sptr); void run_tests(); + + //! Set methods that control which tests are run. + //@{ + void set_do_test_gradient(bool d); + void set_do_test_Hessian_convexity(bool d); + void set_do_test_Hessian_against_numerical(bool d); + //@} + protected: char const * density_filename; shared_ptr > objective_function_sptr; @@ -79,6 +90,47 @@ class GeneralisedPriorTests : public RunTests void run_tests_for_objective_function(const std::string& test_name, GeneralisedPrior& objective_function, shared_ptr target_sptr); + + //! Tests the prior's gradient by comparing to the numerical gradient computed using perturbation response. + void test_gradient(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr); + + //! Test various configurations of the Hessian of the prior via accumulate_Hessian_times_input() for convexity + /*! + Tests the convexity condition: + \f[ x^T \cdot H_{\lambda}x >= 0 \f] + for all non-negative \c x and non-zero \c \lambda (Relative Difference Prior conditions). + This function constructs an array of configurations to test this condition and calls + \c test_Hessian_convexity_configuration(). + */ + void test_Hessian_convexity(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr); + + //! Tests the compute_Hessian method implemented into convex priors + /*! Performs a perturbation response using compute_gradient to determine if the compute_Hessian (for a single densel) + is within tolerance. + */ + void test_Hessian_against_numerical(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr); + +private: + //! Hessian test for a particular configuration of the Hessian concave condition + bool test_Hessian_convexity_configuration(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr, + const float beta, + const float input_multiplication, const float input_addition, + const float current_image_multiplication, const float current_image_addition); + + //! Variables to control which tests are run, see the set methods + //@{ + bool do_test_gradient = false; + bool do_test_Hessian_convexity = false; + bool do_test_Hessian_against_numerical = false; + //@} }; GeneralisedPriorTests:: @@ -86,6 +138,27 @@ GeneralisedPriorTests(char const * const density_filename) : density_filename(density_filename) {} +void +GeneralisedPriorTests:: +set_do_test_gradient(const bool d) +{ + do_test_gradient = d; +} + +void +GeneralisedPriorTests:: +set_do_test_Hessian_convexity(const bool d) +{ + do_test_Hessian_convexity = d; +} + +void +GeneralisedPriorTests:: +set_do_test_Hessian_against_numerical(const bool d) +{ + do_test_Hessian_against_numerical = d; +} + void GeneralisedPriorTests:: run_tests_for_objective_function(const std::string& test_name, @@ -96,46 +169,258 @@ run_tests_for_objective_function(const std::string& test_name, if (!check(objective_function.set_up(target_sptr)==Succeeded::yes, "set-up of objective function")) return; + if (do_test_gradient) + { + std::cerr << "----- test " << test_name << " --> Gradient\n"; + test_gradient(test_name, objective_function, target_sptr); + } + + if (do_test_Hessian_convexity) + { + std::cerr << "----- test " << test_name << " --> Hessian-vector product for convexity\n"; + test_Hessian_convexity(test_name, objective_function, target_sptr); + } + + if (do_test_Hessian_against_numerical) + { + std::cerr << "----- test " << test_name << " --> Hessian against numerical\n"; + test_Hessian_against_numerical(test_name, objective_function, target_sptr); + } +} + + +void +GeneralisedPriorTests:: +test_gradient(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr) +{ // setup images target_type& target(*target_sptr); shared_ptr gradient_sptr(target.get_empty_copy()); shared_ptr gradient_2_sptr(target.get_empty_copy()); info("Computing gradient",3); + const int verbosity_default = Verbosity::get(); + Verbosity::set(0); objective_function.compute_gradient(*gradient_sptr, target); + Verbosity::set(verbosity_default); this->set_tolerance(std::max(fabs(double(gradient_sptr->find_min())), fabs(double(gradient_sptr->find_max()))) /1000); info("Computing objective function at target",3); const double value_at_target = objective_function.compute_value(target); target_type::full_iterator target_iter=target.begin_all(); target_type::full_iterator gradient_iter=gradient_sptr->begin_all(); - target_type::full_iterator gradient_2_iter=gradient_2_sptr->begin_all(); + target_type::full_iterator gradient_2_iter=gradient_2_sptr->begin_all(); // setup perturbation response const float eps = 1e-3F; bool testOK = true; info("Computing gradient of objective function by numerical differences (this will take a while)",3); while(target_iter!=target.end_all())// && testOK) + { + const float org_image_value = *target_iter; + *target_iter += eps; // perturb current voxel + const double value_at_inc = objective_function.compute_value(target); + *target_iter = org_image_value; // restore + const float ngradient_at_iter = static_cast((value_at_inc - value_at_target)/eps); + *gradient_2_iter = ngradient_at_iter; + testOK = testOK && this->check_if_equal(ngradient_at_iter, *gradient_iter, "gradient"); + //for (int i=0; i<5 && target_iter!=target.end_all(); ++i) { - const float org_image_value = *target_iter; - *target_iter += eps; // perturb current voxel - const double value_at_inc = objective_function.compute_value(target); - *target_iter = org_image_value; // restore - const float ngradient_at_iter = static_cast((value_at_inc - value_at_target)/eps); - *gradient_2_iter = ngradient_at_iter; - testOK = testOK && this->check_if_equal(ngradient_at_iter, *gradient_iter, "gradient"); - //for (int i=0; i<5 && target_iter!=target.end_all(); ++i) - { - ++gradient_2_iter; ++target_iter; ++ gradient_iter; - } + ++gradient_2_iter; ++target_iter; ++ gradient_iter; } + } if (!testOK) + { + std::cerr << "Numerical gradient test failed with for " + test_name + " prior\n"; + info("Writing diagnostic files gradient" + test_name + ".hv, numerical_gradient" + test_name + ".hv"); + write_to_file("gradient" + test_name + ".hv", *gradient_sptr); + write_to_file("numerical_gradient" + test_name + ".hv", *gradient_2_sptr); + } +} + +void +GeneralisedPriorTests:: +test_Hessian_convexity(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr) +{ + if (!objective_function.get_is_convex()) + return; + /// Construct configurations + float beta_array[] = {0.01, 1, 100}; // Penalty strength should only affect scale + // Modifications to the input image + float input_multiplication_array[] = {-100, -1, 0.01, 1, 100}; // Test negative, small and large values + float input_addition_array[] = {-10, -1, -0.5, 0.0, 1, 10}; + // Modifications to the current image (Hessian computation) + float current_image_multiplication_array[] = {0.01, 1, 100}; + float current_image_addition_array[] = {0.0, 0.5, 1, 10}; // RDP has constraint that current_image is non-negative + + bool testOK = true; + float initial_beta = objective_function.get_penalisation_factor(); + for (float beta : beta_array) + for (float input_multiplication : input_multiplication_array) + for (float input_addition : input_addition_array) + for (float current_image_multiplication : current_image_multiplication_array) + for (float current_image_addition : current_image_addition_array) { + if (testOK) // only compute configuration if testOK from previous tests + testOK = test_Hessian_convexity_configuration(test_name, objective_function, target_sptr, + beta, + input_multiplication, input_addition, + current_image_multiplication, current_image_addition); + } + /// Reset beta to original value + objective_function.set_penalisation_factor(initial_beta); +} + +bool +GeneralisedPriorTests:: +test_Hessian_convexity_configuration(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr, + const float beta, + const float input_multiplication, const float input_addition, + const float current_image_multiplication, const float current_image_addition) +{ + /// setup targets + target_type& target(*target_sptr); + shared_ptr output(target.get_empty_copy()); + shared_ptr current_image(target.get_empty_copy()); + shared_ptr input(target.get_empty_copy()); + objective_function.set_penalisation_factor(beta); + { + /// Construct an current_image & input with various values that are variations of the target + target_type::full_iterator input_iter = input->begin_all(); + target_type::full_iterator current_image_iter = current_image->begin_all(); + target_type::full_iterator target_iter = target.begin_all(); + while (input_iter != input->end_all()) { - info("Writing diagnostic files gradient" + test_name + ".hv, numerical_gradient" + test_name + ".hv"); - write_to_file("gradient" + test_name + ".hv", *gradient_sptr); - write_to_file("numerical_gradient" + test_name + ".hv", *gradient_2_sptr); + *input_iter = input_multiplication * *target_iter + input_addition; + *current_image_iter = current_image_multiplication * *target_iter + current_image_addition; + ++input_iter; ++target_iter; ++current_image_iter; } + } + + /// Compute H x + objective_function.accumulate_Hessian_times_input(*output, *current_image, *input); + + /// Compute x \cdot (H x) + float my_sum = 0.0; + my_sum = std::inner_product(input->begin_all(), input->end_all(), output->begin_all(), my_sum); + + // test for a CONVEX function + if (this->check_if_less(0, my_sum)) { +// info("PASS: Computation of x^T H x = " + std::to_string(my_sum) + " > 0 and is therefore convex"); + return true; + } else { + // print to console the FAILED configuration + info("FAIL: Computation of x^T H x = " + std::to_string(my_sum) + " < 0 and is therefore NOT convex" + + "\ntest_name=" + test_name + + "\nbeta=" + std::to_string(beta) + + "\ninput_multiplication=" + std::to_string(input_multiplication) + + "\ninput_addition=" + std::to_string(input_addition) + + "\ncurrent_image_multiplication=" + std::to_string(current_image_multiplication) + + "\ncurrent_image_addition=" + std::to_string(current_image_addition) + + "\n >input image max=" + std::to_string(input->find_max()) + + "\n >input image min=" + std::to_string(input->find_min()) + + "\n >target image max=" + std::to_string(target.find_max()) + + "\n >target image min=" + std::to_string(target.find_min())); + return false; + } +} + + +void +GeneralisedPriorTests:: +test_Hessian_against_numerical(const std::string &test_name, + GeneralisedPrior &objective_function, + shared_ptr target_sptr) +{ + if (!objective_function.get_is_convex()) + return; + /// Setup + const float eps = 1e-3F; + bool testOK = true; + const int verbosity_default = Verbosity::get(); + + // setup images + target_type& input(*target_sptr->get_empty_copy()); + input += *target_sptr; // make input have same values as target_sptr + shared_ptr gradient_sptr(target_sptr->get_empty_copy()); + shared_ptr pert_grad_and_numerical_Hessian_sptr(target_sptr->get_empty_copy()); + shared_ptr Hessian_sptr(target_sptr->get_empty_copy()); + + Verbosity::set(0); + objective_function.compute_gradient(*gradient_sptr, input); + Verbosity::set(verbosity_default); +// this->set_tolerance(std::max(fabs(double(gradient_sptr->find_min())), fabs(double(gradient_sptr->find_max()))) /10); + + // Setup coordinates (z,y,x) for perturbation test (Hessian will also be computed w.r.t this voxel, j) + BasicCoordinate<3, int> perturbation_coords; + + // Get min/max indices + const int min_z = input.get_min_index(); + const int max_z = input.get_max_index(); + const int min_y = input[min_z].get_min_index(); + const int max_y = input[min_z].get_max_index(); + const int min_x = input[min_z][min_y].get_min_index(); + const int max_x = input[min_z][min_y].get_max_index(); + + // Loop over each voxel j in the input and check perturbation response. + for (int z=min_z;z<= max_z;z++) + for (int y=min_y;y<= max_y;y++) + for (int x=min_x;x<= max_x;x++) + if (testOK) + { + perturbation_coords[1] = z; perturbation_coords[2] = y; perturbation_coords[3] = x; + + // Compute H(x)_j (row of the Hessian at the jth voxel) + objective_function.compute_Hessian(*Hessian_sptr, perturbation_coords, input); + this->set_tolerance(std::max(fabs(double(Hessian_sptr->find_min())), fabs(double(Hessian_sptr->find_max()))) /500); + + // Compute g(x + eps) + Verbosity::set(0); + // Perturb target at jth voxel, compute perturbed gradient, and reset voxel to original value + float perturbed_voxels_original_value = input[perturbation_coords[1]][perturbation_coords[2]][perturbation_coords[3]]; + input[perturbation_coords[1]][perturbation_coords[2]][perturbation_coords[3]] += eps; + objective_function.compute_gradient(*pert_grad_and_numerical_Hessian_sptr, input); + input[perturbation_coords[1]][perturbation_coords[2]][perturbation_coords[3]] = perturbed_voxels_original_value; + + // Now compute the numerical-Hessian = (g(x+eps) - g(x))/eps + *pert_grad_and_numerical_Hessian_sptr -= *gradient_sptr; + *pert_grad_and_numerical_Hessian_sptr /= eps; + + Verbosity::set(verbosity_default); + // Test if pert_grad_and_numerical_Hessian_sptr is all zeros. + // This can happen if the eps is too small. This is a quick test that allows for easier debugging. + if (pert_grad_and_numerical_Hessian_sptr->sum_positive() == 0.0 && Hessian_sptr->sum_positive() > 0.0) + { + this->everything_ok = false; + testOK = false; + info("test_Hessian_against_numerical: failed because all values are 0 in numerical Hessian"); + } + + // Loop over each of the voxels and compare the numerical-Hessian with Hessian + target_type::full_iterator numerical_Hessian_iter = pert_grad_and_numerical_Hessian_sptr->begin_all(); + target_type::full_iterator Hessian_iter = Hessian_sptr->begin_all(); + while(numerical_Hessian_iter != pert_grad_and_numerical_Hessian_sptr->end_all() && testOK) + { + testOK = testOK && this->check_if_equal(*Hessian_iter, *numerical_Hessian_iter, "Hessian"); + ++numerical_Hessian_iter; ++ Hessian_iter; + } + + if (!testOK) + { + // Output volumes for debug + std::cerr << "Numerical-Hessian test failed with for " + test_name + " prior\n"; + info("Writing diagnostic files `Hessian_" + test_name + ".hv` and `numerical_Hessian_" + test_name + ".hv`"); + write_to_file("Hessian_" + test_name + ".hv", *Hessian_sptr); + write_to_file("numerical_Hessian_" + test_name + ".hv", *pert_grad_and_numerical_Hessian_sptr); + write_to_file("input_" + test_name + ".hv", input); + } + } } void @@ -179,30 +464,52 @@ run_tests() shared_ptr density_sptr; construct_input_data(density_sptr); - std::cerr << "Tests for QuadraticPrior\n"; + std::cerr << "\n\nTests for QuadraticPrior\n"; { QuadraticPrior objective_function(false, 1.F); + set_do_test_gradient(true); + set_do_test_Hessian_convexity(true); + set_do_test_Hessian_against_numerical(true); this->run_tests_for_objective_function("Quadratic_no_kappa", objective_function, density_sptr); } - std::cerr << "Tests for Relative Difference Prior\n"; + std::cerr << "\n\nTests for Relative Difference Prior with epsilon = 0\n"; { - // gamma is default and epsilon is off + // gamma is default and epsilon is 0.0 RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.F); - this->run_tests_for_objective_function("RDP_no_kappa", objective_function, density_sptr); + set_do_test_gradient(true); + set_do_test_Hessian_convexity(true); + set_do_test_Hessian_against_numerical(false); // RDP, with epsilon = 0.0, will fail the numerical Hessian test + this->run_tests_for_objective_function("RDP_no_kappa_no_eps", objective_function, density_sptr); + } + std::cerr << "\n\nTests for Relative Difference Prior with epsilon = 0.1\n"; + { + // gamma is default and epsilon is "small" + RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.1F); + set_do_test_gradient(true); + set_do_test_Hessian_convexity(true); + set_do_test_Hessian_against_numerical(true); // With a large enough epsilon the RDP numerical test will pass + this->run_tests_for_objective_function("RDP_no_kappa_with_eps", objective_function, density_sptr); + } + + std::cerr << "\n\nTests for PLSPrior\n"; + { + PLSPrior objective_function(false, 1.F); + shared_ptr > anatomical_image_sptr(density_sptr->get_empty_copy()); + anatomical_image_sptr->fill(1.F); + objective_function.set_anatomical_image_sptr(anatomical_image_sptr); + // Disabled PLS due to known issue + set_do_test_gradient(false); + set_do_test_Hessian_convexity(false); + set_do_test_Hessian_against_numerical(false); + this->run_tests_for_objective_function("PLS_no_kappa_flat_anatomical", objective_function, density_sptr); } - // Disabled PLS due to known issue -// std::cerr << "Tests for PLSPrior\n"; -// { -// PLSPrior objective_function(false, 1.F); -// shared_ptr > anatomical_image_sptr(density_sptr->get_empty_copy()); -// anatomical_image_sptr->fill(1.F); -// objective_function.set_anatomical_image_sptr(anatomical_image_sptr); -// this->run_tests_for_objective_function("PLS_no_kappa_flat_anatomical", objective_function, density_sptr); -// } - std::cerr << "Tests for Logcosh Prior\n"; + std::cerr << "\n\nTests for Logcosh Prior\n"; { // scalar is off LogcoshPrior objective_function(false, 1.F, 1.F); + set_do_test_gradient(true); + set_do_test_Hessian_convexity(true); + set_do_test_Hessian_against_numerical(true); this->run_tests_for_objective_function("Logcosh_no_kappa", objective_function, density_sptr); } }