Removing requirements for compute_sub_gradient to use sensitivity image

[robbietuk]

This is a fix for #873

Introduces the method

actual_compute_sub_gradient_without_penalty(TargetT& gradient, const TargetT &current_estimate, const int subset_num, const bool do_subtraction)

where the bool do_subtraction (name needs to change) indicates to distributable and the gradient computation whether to do the gradient computation as:
do_subtraction = true (e.g. compute_sub_gradient_without_penalty)

backproj_m [ y/ybar -1 ] 

or do_subtraction = false (e.g. compute_sub_gradient_without_penalty_plus_sensitivity)

backproj_m [ y/ybar ]

[robbietuk]

Sorry, I thought a Draft PR didnt run Travis. It appears that using override results in ALOT of warnings along the lines of:

STIR/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.h:80:5: warning: 'construct_target_ptr' overrides a member function but is not marked 'override' [-Winconsistent-missing-override]

[KrisThielemans]

ok. we'll fix all the override issues in one go as in #860. Let's remove them here then.

[KrisThielemans]

or do_subtraction = false (e.g. compute_sub_gradient_without_penalty_plus_sensitivity)

backproj_m [ y/ybar ] - backproj [ 1 ]

that'll have to be

 backproj_m [ y/ybar ]

[robbietuk]

Corrected the above and removed overrides.

ok. we'll fix all the override issues in one go as in #860. Let's remove them here then.

I think you meant #868

Opinion on the do_subtraction name and/or the actual_compute_sub_gradient_without_penalty() method name?

[robbietuk]

Need to fix the documentation regarding sensitivities:

STIR/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.h

Lines 69 to 72 in 42240de

	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().

[robbietuk]

Interestingly the test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData fails because the numerical gradient does not match the gradient.

I think this is due to a normalisation issue with the new implementation that is caused by STIR's relucatnce to apply normalisation during the gradient computation and previous handling using the sensitivity subtraction.

Let me try and explain. The gradient computation does not want/or need the normalisation (N):

STIR/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx

Lines 1207 to 1217 in dc6f356

	void RPC_process_related_viewgrams_gradient(
	const shared_ptr<ForwardProjectorByBin>& forward_projector_sptr,
	const shared_ptr<BackProjectorByBin>& back_projector_sptr,
	RelatedViewgrams<float>* measured_viewgrams_ptr,
	int& count, int& count2, double* log_likelihood_ptr /* = NULL */,
	const RelatedViewgrams<float>* additive_binwise_correction_ptr,
	const RelatedViewgrams<float>* 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");

This is because it utilised the mathmatical shortcut:
Standard loglikelihood grad = A^T [ Ny/ (NA(x) + b) - 1 ]
STIR master grad = A^T [ y/ (A(x) + b/N) ] - A^T [ 1/N ]

A is the projector operator and b is an additive term (with the state of normalisation unknown, it isnt that important for this). This allows for the precomputation of the constant terms: additive (b/N) and sensitivity (A^T [ 1/N ]). Then, the gradient computation was simple and the RPC_process_related_viewgrams_gradient does not apply the normalisation to either the expected or the measured data.

If this is the case, this PR modification might need some tweaking. Rather than subtracting ones, it should subtract 1/N. Correct?

[robbietuk]

Test passes on my end. I get the same result (ish)

gradient:

difference between compute_sub_gradient_without_penalty_plus_sensitivity() - subset_sensitivity and gradient.hv

Edit: This was fixed by fa9e87e

[robbietuk]

Test failing on OSSPS in recon_test_pack

STIR/recon_test_pack/run_tests.sh

Lines 231 to 246 in 42240de

	echo
	echo -------- Running OSSPS with a quadratic prior --------
	echo Running ${INSTALL_DIR}OSSPS
	${MPIRUN} ${INSTALL_DIR}OSSPS OSSPS_test_PM_QP.par 1> OSSPS_PM_QP.log 2> OSSPS_PM_QP_stderr.log
	echo '---- Comparing output of OSSPS subiter 8 (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_8.hv -a
	${INSTALL_DIR}compare_image -r 1 test_image_OSSPS_PM_QP_8.hv my_test_image_OSSPS_PM_QP_8.hv
	then
	echo ---- This test seems to be ok !;
	else
	echo There were problems here!;
	ThereWereErrors=1;
	fi

This PR

The error output

-------- Running OSSPS with a quadratic prior --------
Running OSSPS
---- Comparing output of OSSPS subiter 8 (should be identical up to tolerance)
Running compare_image

Maximum absolute error = 0.00227004
Maximum in (1st - 2nd) = 0.00227004
Minimum in (1st - 2nd) = -0.000148283
Error relative to sup-norm of first image = 1.01647 %

Image arrays deemed different
(tolerance used: 0.2 %)


Maximum absolute error = 0.00131911
Maximum in (1st - 2nd) = 0.00131911
Minimum in (1st - 2nd) = -1.35558e-05
Error relative to sup-norm of first image = 0.590663 %

Image arrays deemed different
(tolerance used: 0.05 %)

There were problems here!

The OSSPS reconstruction is fine, just the reconstructed image and reference are ~1% different.

STIR Master

The difference between the reference in STIR master are of magnetude e-04 % / e-05 %.

Discussion

I think this is because the OSSPS reconstuction parameter file using use subset sensitivities:=0.
In STIR master, :
this->objective_function_sptr->compute_sub_gradient(*numerator_ptr, current_image_estimate, subset_num); would be subtracting the subset sensitvity from backproj[ y/ybar] but the subset sensitivity is the total_sensitivity/ num_subset.
In this PR, the gradient will subtract the actual subset sensitivity (or the equivalent.
Therefore, I think the issue is the test and the reference image (test_image_OSSPS_PM_QP_8.hv) is no longer correct for this objective function gradient in the OSSPS reconstruction.

Difference between master image and reference is virtually 0
Difference between PR image and reference:

Most of the error is in the low sensitivity region, which definitely makes sense if the total sensitivity was subtracted.

Proposed Solution

Replace the OSSPS reference image.

[robbietuk]

Now that I have confirmed that there is a difference between OSSPS when use subset sensitivities:= 0 or 1, the difference between the image reconstructed using STIR master with use subset sensitivities:=1 and this PR is:

This image is coronal view (top/bottom are the end rings). Clearly there is an end plane issue but the majority of the interior is good.

I tested to see if this was a subset sensitivity issue in STIR master but comparing the resulting images output by the test when use subset sensitivities:= 0 and 1
There are minor issues but the lack of end plane difference suggests that it is a result of this PR.

This is affected by changing from zero end planes of segment 0:= 1 to zero end planes of segment 0:= 0 but it does not resolve the issue, as shown here.

Maybe this is a normalisation subtraction issue in the PR?

The difference is subtle and does occur over 8 OSSPS sub-iterations with 4 subsets

[robbietuk]

Major changes since last comments:

• Resolved the endplane issue. It was a normalisation issue.
• Updated the OSSPS recon_test_pack reference image to this PR (previously failed because the reference was computed with use_subset_sensitivities =0 which resulted in the wrong objective function.
• Added OSSPS recon_test_pack test to compare the output using subset senstitivites and not. These should be identical (expect for list mode right now [not tested in recon_test_pack], which I am not sure how to fix).
• Added compute_sub_gradient_without_penalty documentation

Marking as ready for review (pending travis). There are still a few issues:

• List mode compute_sub_gradient_without_penalty() has been made to error when use_subset_sensitivities=true because this is the only call now that does not use the "standard" objective function. Keep the error or let it run normally? Maybe add lots of warning? up to you @KrisThielemans
• In distributable_compute_gradient() the normalisation_sptr should be parsed in both cases. This is probably safer to do anyway.

Edit: Travis passed on previous major commit

[robbietuk]

Ready for review