Improvements to scatter estimation interpolation

[markus-jehl]

In scatter estimation the projection data are interpolated during upsample_and_fit_scatter_estimate, which calls interpolate_projdata. This function uses sample_function_on_regular_grid for Cylindrical scanners and interpolate_axial_position for BlocksOnCylindrical scanners. In the axial direction, this results in very unsmooth scatter estimates for BlocksOnCylindrical scanners:

Therefore, it would be nice to rewrite the interpolate_projdata function to take the scanner geometry into account when calling the interpolator (basically replacing sample_function_on_regular_grid with e.g. sample_function_on_ProjData.

[markus-jehl]

@NikEfth : @KrisThielemans mentioned that you are currently working on something related, so it would be good to hear your take on this!

[NikEfth]

Hi @markus-jehl and @KrisThielemans,

In the past, I was working on energy-based scatter estimation. The method returned down-sampled 3D sinograms which I then had to interpolate back to the original size.

I had tried many things but the one that worked the best was to split the interpolation into 2 steps. First, I interpolate the views and tangs from the downsampled to the full geometry in an in-segment fashion (not to be confused with direct, as obliques were included, too). Then I convert the partially upsampled sinograms to michelograms and interpolate between rings. The results are good and you are welcome to give it a try.
The only thing you would need to do is to change the flag to do a 3D scatter simulation.

Of note is that in the 2D the typical output is the scatter and additive sinograms. Here the scatter is not exported as essentially the scaled are scattered. So to regenerate your figure above use the scaled from the folder extras.

I could make this more consistent tbh.

[NikEfth]

I am sorry, In my previous reply, I forgot to mention the BlocksOnCylindrical case.
So, right now, I left that case as is, because was not sure about @KrisThielemans plans.
You can try to interpolate them quite easily, again by activating the 3D scatter estimation, and commenting out the line in interpolate_3d that checks for BlocksOnCylindrical. Then it should move forward with the interpolation.

I am looking forward to your results!

[markus-jehl]

Hi Nikos,
Thanks for the quick reply! Testing it now (not sure how long this will run for :-) ) - but at least that gives me time to go through the code. Will let you know once I have results.

[markus-jehl]

Is this the branch I should be using? https://github.com/NikEfth/STIR/tree/scatter_3d

Have been pulled away on something else and will try this in spare moments - expect a bit of a delay until I have results :-)

[KrisThielemans]

This will be #1156 presumably.

As far as I can see, this PR doesn't resolve the issue of gaps in axial direction as it still uses

STIR/src/buildblock/interpolate_projdata.cxx

Line 444 in 9ed30a6

sample_function_on_regular_grid(sino_2D_out, proj_data_interpolator, offset, step);

Of course, the NeuroLF doesn' t have any, so it might be sufficient for @markus-jehl

[markus-jehl]

@KrisThielemans: I've started looking into interpolate_projdata.cxx a bit and trying to work out the best way to implement this. The approach in the description of this issue would certainly be the nicest, but that only works if there is an interpolator that can work on irregular input data (currently it uses BSplinesRegularGrid) - is there one and how well do you think it will work if there are only around 7 slices in the input proj data?

Alternatively, I can just add some axial smoothing to interpolate_axial_position.cxx, which seems fairly straight-forward.

Also, should this line not be adding the offset instead of subtracting?

STIR/src/include/stir/numerics/sampling_functions.inl

Line 34 in bf0366a

relative_positions[1]= index_out[1] * step[1] - offset[1] ;

In all other places in this function, the offset is added.

@NikEfth: I tried your branch and setting 3D, but all my RAM was being filled up and it took forever. Probably I would have needed to reduce the segment numbers etc. to get it working, but since all I need is a bit of axial smoothing on the direct scatter estimate I thought it's easier to do that directly.

[NikEfth]

Sigh. That's a pity I was really looking forward to these results.

[KrisThielemans]

The approach in the description of this issue would certainly be the nicest, but that only works if there is an interpolator that can work on irregular input data

but the way we construct the downsampled scanner is "regular" (i.e. no gaps there), so I still think this should work. Of course, this would be inappropriate for scanners with huge gaps (as then we'd potentially be computing scatter inside a large gaps, and never use it again).

In any case, the (spatially invariant) smoothing approach assumes that the output sampling is regular, which is fine for you, but not for other scanners.

should this line not be adding the offset instead of subtracting?

I will have to check that later (but I hope this stuff is tested elsewhere)

[NikEfth]

this would be inappropriate for scanners with huge gaps (as then we'd potentially be computing scatter inside a large gaps, and never use it again).

This is a problem I had to face when a small number of blocks are off, too.
So, I create an extra sinogram with the 0 in the dead blocks and gaps and multiplied with the upsampled before tail fitting. To get the correct number of events.
I have this code lying somewhere, too.

[KrisThielemans]

I have this code lying somewhere, too.

as far as I can see, the only thing we need is changing the interpolator to cope with non-uniform output. All the rest will still work.

[markus-jehl]

I have this code lying somewhere, too.

as far as I can see, the only thing we need is changing the interpolator to cope with non-uniform output. All the rest will still work.

on it :-)