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

Issue/1161 improvements to scatter estimation interpolation #1172

Conversation

markus-jehl
Copy link
Contributor

@markus-jehl markus-jehl commented Mar 7, 2023

Linked to issue: #1161

  • Fixed InverseSSRB for BlocksOnCylindrical by completely rewriting it.
  • Modified interpolate_projdata significantly, to not only perform smoothing for Cylindrical but also BlocksOnCylindrical (in axial direction)
  • Removed reliance on extend_segment_in_views, which (if at all) only worked for Cylindrical. Instead, extending it very simply in the code directly. Also, previously some dimensions were not extended while others were - all this is tidied up now.
  • Still to do: add test

@markus-jehl
Copy link
Contributor Author

Scatter estimation from simulated Moliner phantom data before and after, plus match with ground truth:
image
image
image

@markus-jehl
Copy link
Contributor Author

BlocksOnCylindrical, quadratic b-spline:
image
BlocksOnCylindrical, linear b-spline:
image
BlocksOnCylindrical: artefact when not extending proj data for interpolator
image

@markus-jehl
Copy link
Contributor Author

Cylindrical, quadratic b-spline:
image
Cylindrical, linear b-spline:
image

@markus-jehl
Copy link
Contributor Author

Cylindrical and BlocksOnCylindrical InverseSSRB result:
image
image

Markus Jehl added 3 commits March 7, 2023 16:12
@KrisThielemans
Copy link
Collaborator

AppvVeyor complains with

C:\projects\stir\src\buildblock\extend_projdata.cxx(176,130): error C2065: 'M_PI': undeclared identifier 

Will also need an update in the release notes

@markus-jehl
Copy link
Contributor Author

Also fixes #1177

@markus-jehl
Copy link
Contributor Author

Also fixes #1178

…check to ensure views and tangential bins are the same number for BlocksOnCylindrical.
Copy link
Collaborator

@KrisThielemans KrisThielemans left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't reviewed interpolate_projdata.cxx very well yet. It looks reasonable but due to white-space changes, I'm not sure if it does the same as before.

src/test/test_interpolate_scatter.cxx Outdated Show resolved Hide resolved
src/test/test_interpolate_scatter.cxx Outdated Show resolved Hide resolved
Comment on lines 179 to 180
if (abs(segment.get_proj_data_info_sptr()->get_phi(Bin(0, segment.get_proj_data_info_sptr()->get_max_view_num(), 0, 0)) -
segment.get_proj_data_info_sptr()->get_phi(Bin(0, segment.get_proj_data_info_sptr()->get_min_view_num(), 0, 0))) > _PI)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This doesn't really test for 360, but larger than 180. To cope with SPECT, I'd prefer to make this test more accurate

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is actually not trivial, because for BlocksOnCylindrical scanners phi doesn't change continuously, so I can't just compare the first and last phi to _PI. Can you think of a check for specifically 180° and 360° that works for both BOC and Cylindrical?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

for BlocksOnCylindrical scanners phi doesn't change continuously

I guess you mean it isn't regularly sampled. True. Would something like this work?

const float phi_diff = max_phi - min_phi;
const float average_phi_sampling = phi_diff/(get_num_views()-1);
const float approx_phi_range = average_phi_sampling * get_num_views();
if (abs(phi_range - 2*_PI) < average_phi_sampling)
...

a bit elaborate of course

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

image
unfortunately this doesn't work for BlocksOnCylindrical

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ah well. a factor ~4 off. I guess this will depend on the size of the gaps etc.

It brings us to an interesting question. The output of extend_segment is used as input for interpolation. If the views are very different, then you could argue that we shouldn't be using the "wrap-around" ones anyway.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is true, if the views are very different it would make more sense to just extend with the current view. However, for symmetric BlocksOnCylindrical scanners it would make sense to wrap around, since this would be exactly the same view difference that we also interpolate in other places of the sinograms.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have now set it to use 5* average_sampling as the cutoff, and to fill with nearest neighbour and print a warning if it is further away from either 180° or 360°

src/buildblock/extend_projdata.cxx Outdated Show resolved Hide resolved
src/buildblock/inverse_SSRB.cxx Show resolved Hide resolved
src/buildblock/inverse_SSRB.cxx Outdated Show resolved Hide resolved
src/include/stir/extend_projdata.h Show resolved Hide resolved
src/include/stir/extend_projdata.h Outdated Show resolved Hide resolved
src/buildblock/interpolate_projdata.cxx Outdated Show resolved Hide resolved
src/buildblock/interpolate_projdata.cxx Show resolved Hide resolved
@KrisThielemans
Copy link
Collaborator

By the way, is interpolate_axial_position now obsolete? If so, please deprecate

@KrisThielemans
Copy link
Collaborator

Let's remove the obsolete warning at https://github.com/markus-jehl/STIR/blob/d9aa91bcd260efcf57152a40becd4770eff4a148/src/buildblock/interpolate_projdata.cxx#L190-L191

@KrisThielemans
Copy link
Collaborator

@markus-jehl
Copy link
Contributor Author

Let's remove the obsolete warning at https://github.com/markus-jehl/STIR/blob/d9aa91bcd260efcf57152a40becd4770eff4a148/src/buildblock/interpolate_projdata.cxx#L190-L191

actually, the relevant code in Cylindrical looks rather suspicious https://github.com/markus-jehl/STIR/blob/d9aa91bcd260efcf57152a40becd4770eff4a148/src/buildblock/interpolate_projdata.cxx#L236-L239 Did we ever test this?

Yes, the tests I added do interpolation of Cylindrical data with view offsets and everything looks as expected. What part of the code looks suspicious to you?

@KrisThielemans
Copy link
Collaborator

thanks @markus-jehl can you also remove the obsolete warning? I would need stating in the release notes as well.

@KrisThielemans
Copy link
Collaborator

@danieldeidda @markus-jehl given the back-and-forth a bit on this PR, ok to squash-merge?

@danieldeidda
Copy link
Collaborator

Ok for me

@markus-jehl
Copy link
Contributor Author

The warning is already removed, but will add this to the release notes.

@KrisThielemans
Copy link
Collaborator

Clean version #1181 merged

@markus-jehl markus-jehl deleted the issue/1161-improvements-to-scatter-estimation-interpolation branch April 11, 2023 07:21
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

Successfully merging this pull request may close these issues.

None yet

3 participants