BlocksOnCylindrical get_sampling_in_m returns ring spacing instead of taking gaps into account
markus-jehl opened this issue · 6 comments
As part of the scatter speedup work with downsampled scanners (#1291), I noticed that there was always a large artifact at the bottom of the scatter estimate (~15% error), irrespective of how much downsampling was used. This turned out to be caused by the interpolation step that is done in the final step of the scatter estimation (that interpolation is actually not required, since the scatter is already upsampled to the original scanner dimensions). There, the input and output are direct sinograms of the same size and the check inside upsample_and_fit_scatter_estimate, whether the input is regularly sampled succeeds because get_sampling_in_m trivially returns the ring spacing. However, since here the input is the full BlocksOnCylindrical scanner and no longer the downsampled scanner, it should actually fail because there are gaps between axial blocks:
As a result, the bottom of the simulated cylinder is now ~1mm off after this second interpolation step and causes the horizontal artifact:
@KrisThielemans @danieldeidda
This is the problem:
STIR/src/include/stir/ProjDataInfoGeneric.inl
Lines 88 to 94 in 8ced2d7
Compared to what is done for Cylindrical scanners:
STIR/src/buildblock/ProjDataInfo.cxx
Lines 103 to 104 in 8ced2d7
Why are we not using the same implementation for BlocksOnCylindrical? This would have worked:
I see a number of other TODOBLOCK in *Generic.inl ...
Checking the history here, @danieldeidda put those lines in, simplifying @roemicha cde7f2d's copy of the cylindrical case. The simpler code would do exactly the same as the cylindrical case I believe (if span=1). (I guess both cases were implemented in a specific way to speed it up a bit, although not so sure if it actually does).
I see indeed no reason to overload get_sampling_in_m et al. in ProjDataInfoGeneric. Note however that due to our still crazy hierarchy (ProjDataInfoGeneric is still derived from ProjDataInfoCylindrical as opposed to other way around), removing get_sampling_in_m() for Generic won't work. It would have to explictly call the ProjDataInfo version. That would lead some even uglier code, so my current impression is to remove all overloads of get_sampling_in_m/t/theta (I didn't check s and phi)
Another question I suppose is if the down/upsampling code should really call this function. I suppose maybe it needs to for the interpolation, but it certainly shouldn't have caused a drift.
I agree that the cleanest would be to remove all overloads. When the crazy hierarchy is fixed then the implementation can move from the Cylindrical to the Generic projdata. Will give it a try.
Currently the interpolation code is calling this function to check that the input projdata are homogeneously sampled - otherwise the index conversion would be much more complex and computationally demanding.
Currently the interpolation code is calling this function to check that the input projdata are homogeneously sampled - otherwise the index conversion would be much more complex and computationally demanding.
yes. That check could be written in terms of get_m which would be a bit clearer I guess, but fine. However, as in the scatter code the downsampled scanner actually does have equidistant sampling, isn't the problem there that it has the wrong value? (I'm confused)
Currently the interpolation code is calling this function to check that the input projdata are homogeneously sampled - otherwise the index conversion would be much more complex and computationally demanding.
yes. That check could be written in terms of
get_mwhich would be a bit clearer I guess, but fine. However, as in the scatter code the downsampled scanner actually does have equidistant sampling, isn't the problem there that it has the wrong value? (I'm confused)
Yes, the check could also just use get_m and then we wouldn't even have to fix the get_sampling_in_m function. Exactly, in the scatter estimation we have no issues currently. The issue only arises once we complete the final iteration and enter that last section where we call update_and_fit_scatter_estimate another time. This time both input and output ProjData are the same size and both have gaps. As mentioned on #1291, that second time we don't actually need to call interpolate_projdata at all and I fixed it on that PR already by replacing it with inverse_SSRB and normalisation.undo().
ok. sorry.
Obviously best to fix this one in any case. Looks like we can do it in 2 separate PRs then, which is great.