ScatterEstimation convergence problem fix

[markus-jehl]

Fix for #1112.

[KrisThielemans]

looks good. can you check what happens if the initial image isn't initialised? (presumably taken care of by get_initial_data_sptr though)

[markus-jehl]

Just ran it and it works fine, as expected!

[KrisThielemans]

great. It'd be reassuring if you or @emikhay could run some more scatter estimations on your previous data with this branch before we merge. Also, please add something to the release notes (`[ci skip]')

[markus-jehl]

Triple-checked it with simulated data and BPET reconstructions with different normalisations - everything looks good.

[KrisThielemans]

I ran examples/Siemens_mMR/process_data.sh on the corresponding Zenodo data. Below is a screenshot where the bottom image is the difference images (displayed with colour scale [-1%,+1%] of the max used for the top images).

ROI values are

# Stats for Study: Unknown	Generated on: Wed Feb 15 23:15:14 2023
# ROI       	Data Set    	       Frame	Duration (s)	   Midpt (s)	        Gate	Gate Time (s)	      Median	        Mean	         Var	     Std Dev	         Min	         Max	 Size (mm^3)	Frac. Voxels	      Voxels
L           	PR_final_activity_image_42	           0	     500.000	     250.000	           0	       0.000	    0.126365	    0.132625	  0.00215737	   0.0464475	  0.00735384	    0.322484	       18813	     2127.94	        2592
L           	master_final_activity_image_42	           0	     500.000	     250.000	           0	       0.000	    0.126569	    0.132807	   0.0021592	   0.0464672	  0.00739306	    0.322794	       18813	     2127.94	        2592
L           	PR - master 	           0	     500.000	     250.000	           0	       0.000	-0.000174373	-0.000177654	 6.42112e-09	 8.01319e-05	-0.000559345	 9.48533e-05	     18810.9	     2244.50	        2720
B           	PR_final_activity_image_42	           0	     500.000	     250.000	           0	       0.000	   0.0113327	   0.0131966	  6.5529e-05	    0.008095	 0.000809588	   0.0854493	       92693	    10484.47	       12064
B           	master_final_activity_image_42	           0	     500.000	     250.000	           0	       0.000	   0.0111794	   0.0130344	 6.43831e-05	  0.00802391	 0.000799468	   0.0847749	       92693	    10484.47	       12064
B           	PR - master 	           0	     500.000	     250.000	           0	       0.000	 0.000142739	  0.00016161	  7.6954e-09	 8.77234e-05	-5.27129e-07	 0.000925265	       92656	    11055.62	       12558

Showing that this can give a ~1% difference (in the background). It's impossible to know which one is more correct I think, but it is probably a larger difference than people would prefer.

Sadly, I think we should provide an option to be able to do the same thing as before (but leave it off by default).

[markus-jehl]

Hmm... I guess modifying how many iterations the scatter and the underlying OSMAPOSL performs will have a similar impact then? Because all this change is doing is offset the starting point of OSMAPOSL, making it slower to converge. Would it be possible to suggest bumping up the OSMAPOSL iterations as the "legacy mode" you are proposing, instead of actually introducing a new setting?

[KrisThielemans]

it's a pain... but I think user experience isn't great if we release 5.2 which gives scatter estimation differences.

Of course, scatter estimation results depend a lot on specific settings used. The Siemens example uses 21 OSEM updates per scatter iteration, so not so small.

STIR/examples/Siemens-mMR/scatter_and_recon.sh

Lines 80 to 81 in 2e68f4d

	scatter_recon_num_subiterations=21
	scatter_recon_num_subsets=21

Image values (max around 0.12) are also not so far off from 1. So that means people could see bigger differences in their data.

I suggest therefore introducing a new variable and method restart_reconstruction_every_scatter_iteration or similar (defaulting to false). I think it's only the "second" fill(1) that we'd need to optionally disable, as the first one is a bug (we can set the initial image, but was ignored I guess). It'd have to be documented of course. So not so trivial sadly. Feel free to decline doing that!

[robbietuk]

I suggest therefore introducing a new variable and method restart_reconstruction_every_scatter_iteration or similar (defaulting to false).

Something to consider with this methodology is the persistent OSEM image might begin to converge and noisy. I am not sure whether this would happened or what its effect on scatter estimation would be.

[markus-jehl]

You raise a good point: it sounds to me like the total sum of OSEM iterations is the relevant factor. If we reset the image, then we need to perform this total sum each scatter estimation iteration. Otherwise, we can divide the sum by the scatter estimation iterations and only do this fraction in each iteration (saving a lot of time).

[markus-jehl]

@KrisThielemans : happy to add this setting (especially if it is defaulted to false ;-) ).

[KrisThielemans]

I'm fairly sure that we heavily postfilter (necessary for the downsampling anyway), so I don't think noise increase is a big issue. Of course, this is a user-tunable parameter, so it likely could use some more doc.

[robbietuk]

You raise a good point: it sounds to me like the total sum of OSEM iterations is the relevant factor. If we reset the image, then we need to perform this total sum each scatter estimation iteration. Otherwise, we can divide the sum by the scatter estimation iterations and only do this fraction in each iteration (saving a lot of time).

Yes, but this is getting into the specific implementation that is probably best left to the user. Adding restart_reconstruction_every_scatter_iteration (default: false) sounds like a good idea to me and allows a user to configure the estimation to their desire.

[NikEfth]

I am just catching up with this.
Before you accept it could I ask for see some of the reconstructed images from each iteration?
I mean the downsampled.

Are you sure that the scaled normalization factors do actually normalize your data?

The issue with the normalization in scatter estimation is something that still gets under my skin.
So let's say that at some point you have measurements for normalization (these sum to 10^8 or more) and have gaps if cylindrical. This is "A" and you should not use it.

If you invert "A" and put in the gaps 10^12. This is "B" and you should use this.
If you have component-based normalization estimated by STIR from "A" with 0 in the gaps, this is "C" and you should not use this!

If you invert the "C" and put 10^12 in the gaps this is "D" and you should not use this!

Edit (I forgot one case):
If you have "C" and just put 10^12 in the gaps (without inversion) this is "E" and I have used it in past with good results, but I prefer "B", instead.

---

So if in your case, your scaling go from "B" to "A" you should not do that.

I don't think this option is correct.

Edit:
Basically, use only "B".

[KrisThielemans]

@NikEfth, I'll try to answer this, but it might be hard to write this up. In any case, I can't see how this is related with this PR? The problem that is resolved with this PR is if you scale the norm factors (without gaps), the image should scale appropraitely, but the (upsampled) scatter estimate should be the same. With the existing code, it'd only do that if you let OSEM iterate for a long time (as its convergence rate depends on the initial scale of the image as there is a background term).

Of course, still good to see the intermediate images.

[NikEfth]

the image should scale appropraitely,

I agree. But if you scale the norm factors as described in #1112, it might not be the case anymore.

[KrisThielemans]

the image should scale appropraitely,

I agree. But if you scale the norm factors as described in #1112, it might not be the case anymore.

with this PR, it will. Obviously, there's still the issue that without providing an initial image, OSMAPOSL will initialise with 1, which might be a terrible choice for your data, and that it will still take a while to recover the correct scale.

[KrisThielemans]

On your previous comments

So lets say that at some point you have measurements for normalization (these sum to 10^8 or more) and have gaps if cylindrical. This is "A" and you should not use it.

What do you mean with “you should not use it”? How? If you give measured data from a normalisation scan as input to BinNormalisationFromProjData, that’ll indeed be wrong (independent of any gaps) as that expects the “normalisation factors”

If you invert "A" and put in the gaps 10^12. This is "B" and you should use this.

Yes, because we cannot (reliably) store infinite normalisation factors. (Note of course that with this “direct normalisation” approach, noise in the measurements is an issue, but that’s not what we're talking about here)

If you have component-based normalization estimated by STIR from "A" with 0 in the gaps, this is "C" and you should not use this!

What do you mean? How would you use it? If you use apply_normfactors3D to obtain a projdata that you’ll use in a BinNormalisationFromProjData, you have to be careful that it “normalises” a projdata of 1s. I’m hoping it’ll do the correct thing (i.e. put large values in the gaps). (I'm pretty sure it does, as that's what we're using for the mMR).

If you invert the "C" and put 10^12 in the gaps this is "D" and you should not use this!

Obviously by now I’m lost what you mean with C and D.

(By the way, we now have a BinNormalisationPETFromComponents were most component-based functionality is moved. We cannot parse it yet from within a .par file (I'd like to change the output of the ML_norm* stuff first).)

[NikEfth]

What do you mean with “you should not use it”? How? If you give measured data from a normalisation scan as input to BinNormalisationFromProjData, that’ll indeed be wrong (independent of any gaps) as that expects the “normalisation factors”

Well, I have seen people using this file, so I tried to cover all bases.

Yes, because we cannot (reliably) store infinite normalisation factors. (Note of course that with this “direct normalisation” approach, noise in the measurements is an issue, but that’s not what we're talking about here)

Indeed, it is the direct normalization and I recommend using this for the scatter estimation.
Of course, how this was obtained matters. But for the more part, the noise in the downsampled images is not a huge problem.
In addition, the estimation of the component-based sinogram considers that the events are true. Which is not true during the estimation.

What do you mean? How would you use it? If you use apply_normfactors3D to obtain a projdata that you’ll use in a BinNormalisationFromProjData, you have to be careful that it “normalises” a projdata of 1s. I’m hoping it’ll do the correct thing (i.e. put large values in the gaps). (I'm pretty sure it does, as that's what we're using for the mMR).

Indeed, but for true events. So the component-based worked better for me in the final reconstruction.

It is my understanding that this discussion started when @markus-jehl tried to scale the norm factors by a large value. And I believe that this was the problem with the reported behaviour.

[KrisThielemans]

I don't disagree with your statements on use of direct/component-based norm factors and scatter, but I don't think this has anything to do with this PR.

It is my understanding that this discussion started when @markus-jehl tried to scale the norm factors by a large value.

indeed.

And I believe that this was the problem with the reported behaviour.

Not sure what you mean with "this". In any case, this PR resolved that problem.

[KrisThielemans]

I've re-run the Siemens NEMA example data. Unfortunately this has uncovered another problem. For that data, both master and this PR give large oscillations in the estimated image and therefore scatter estimate. A bit more in #1163.

I mention this here as it would influence the difference between this PR and master that I saw above.

[KrisThielemans]

@NikEfth we're going to merge this soon, unless you have any concrete proof that this is buggy?