Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

BlocksOnCylindrical get_sampling_in_m returns ring spacing instead of taking gaps into account #1396

Open
markus-jehl opened this issue Mar 8, 2024 · 6 comments

Comments

@markus-jehl
Copy link
Contributor

markus-jehl commented Mar 8, 2024

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:
image

As a result, the bottom of the simulated cylinder is now ~1mm off after this second interpolation step and causes the horizontal artifact:
image

@markus-jehl
Copy link
Contributor Author

markus-jehl commented Mar 8, 2024

@KrisThielemans @danieldeidda
This is the problem:

ProjDataInfoGeneric::get_sampling_in_m(const Bin& bin) const
{
// TODOBLOCK
return get_scanner_ptr()
->get_ring_spacing(); // TODO currently restricted to span=1 get_num_axial_poss_per_ring_inc(segment_num);
// return get_axial_sampling(bin.segment_num());
}

Compared to what is done for Cylindrical scanners:
return abs(get_m(Bin(bin.segment_num(), bin.view_num(), bin.axial_pos_num() + 1, bin.tangential_pos_num()))
- get_m(Bin(bin.segment_num(), bin.view_num(), bin.axial_pos_num() - 1, bin.tangential_pos_num())))

Why are we not using the same implementation for BlocksOnCylindrical? This would have worked:
image

@KrisThielemans
Copy link
Collaborator

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.

@markus-jehl
Copy link
Contributor Author

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.

@KrisThielemans
Copy link
Collaborator

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)

@markus-jehl
Copy link
Contributor Author

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)

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

@KrisThielemans
Copy link
Collaborator

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants