Flying Focal Spot for FanBeam and ConeBeam geometries#1509
Conversation
|
Checking updated PR...
Comment last updated at 2020-09-05 19:33:41 UTC |
kohr-h
left a comment
kohr-h
left a comment
There was a problem hiding this comment.
I highly appreciate the effort to support FFS, it's a long-standing open issue we had while working with the medical data. As usual, code and documentation are of good quality, as far as I have seen.
Regarding the implementation, I think we need to take a step back and see if we can't come up with something that is more generic than this approach. I feel that the principle is generic ("apply some shifts at each angle") but the implementation is specifically tailored to FFS, so that the next similar thing needs another extension.
What I could imagine is to accept a function like src_shift_func(angle, geometry) that produces the shifts for a given angle. We could supply the flying_focal_spot function as an example (in this module) and point to it. The good thing about this is that it allows other functions that a user can come up with.
Opinions?
| center_to_src_init = center_to_src_init + ffs_shift | ||
| # broadcasting to perform matrix multiplication "manually", | ||
| # since existing numpy functions do cross product along the outer | ||
| # dimensions, which we don't need here |
There was a problem hiding this comment.
I think you can achieve this with numpy.matmul (which I didn't know either when I wrote similar code).
There was a problem hiding this comment.
I think einsum is usually the most beautiful way to do it.
|
Hi! Thank you for the review @kohr-h! The main issue with a function of the type src_shift_func(angle, geometry) is that it is not clear how to define a shift for an arbitrary angle. Basically, FFS is defined for a sequence of angles that correspond to angular partition. One could use interpolation to define it for any angle, but it seems a bit strange to me to do it just for flexibility of the interface. Or maybe it is not strange? What do you think? |
|
That's exactly the main concern, yes. However, you needed to define FFS for arbitrary angle, too, so it's a concern no matter what. |
|
Why is it necessary to define FFS for arbitrary an angle? |
In practice it isn't, but you still had to do it because the methods require it: odl/odl/tomo/geometry/conebeam.py Lines 405 to 407 in 7fef8ad Other kinds of shifts would look similar. |
|
Hi! I have been thinking about this and I am not convinced that FFS should be a function of the rotation angle. Because I imagine that in practical CT scanners both FFS and rotation angle depend on time. Since we don't have a notion of time in the geometry, we could use the rotation angle as a proxy for time, but I don't see how this would make life easier for people who are trying to define arbitrary shifts. Secondly, even if we define FFS(angle), I don't think that we should allow an arbitrary angle. Firstly, I was concerned only about introducing unnecessary complexity in the code. Now I am thinking that we actually don't know how the source moves between the shifts, e.g. if the speed is uniform, if it is not, doing linear interpolation is wrong and can lead to "silent errors". To circumvent this issue, we could define FFS only for angles that belong to discretization and rise "Not implemented" or "Invalid argument error", if the provided angle does not belong to angular partition. |
That's a valid concern. However, note that geometries have a notion of "motion parameter", which usually is, but doesn't have to be, a rotation angle. It can be anything that parametrizes the part of the geometry that is not the detector coordinate. So it could easily be time. The
The problem with this definition of FFS is that it "sticks" to the time steps, i.e., if you decide to double the frequency of the projections being taken, the frequency of the FFS will also double. It's not a function of continuous time, but of discrete time points. If that's the correct behavior, then okay.
That's correct, and a very valid concern.
That suggestion is impractical, though, because you need to compare floating point values to do the check, and floating point numbers are practically never equal. Checking for "almost equal" is also tricky, in particular around 0. So raising errors for invalid input will inevitably generate annoyance. The best thing I can think of here is to do nearest neighbor interpolation, i.e., to always "snap" to the known values. That gives correct results for the given angles, and still allows arbitrary input. |
|
Sorry, it took such a long time, I have added a function flying_focal_spot in a separate file, which just snaps to the known values. It is possible to provide this function as an argument to geometry. |
|
What is the status of this issue? |
@JevgenijaAksjonova, could you provide an update please? |
|
Oh dear. If I remember correctly this is just a final review short of being merged 😬 I'll give it a fresh look asap. |
|
It would be very convenient, it appears that this new Mayo clinic data has it for some of the scans! |
|
Oh did Mayo finally release their data, do you have a link? I was buggering them for years. |
|
Mayo Clinic announced the public availability of 299 patient CT cases containing head, chest and abdomen exams from Siemens and GE scanners. All cases have routine dose and simulated low dose projection data, full dose images, and clinical data identifying all pathology. They appreciate receiving a copy of any publications that make use of these data. Detailed information and a link to the data library are available here: https://ctcicblog.mayo.edu/hubcap/patient-ct-projection-data-library/. |
|
I did manage to read their data in python, but they provide this FFS in some of them. Would be good if this gets implemented, then I can publish some code giving a baseline. |
kohr-h
left a comment
kohr-h
left a comment
There was a problem hiding this comment.
Some suggestions and comments, first round. Don't have more time right now, will come back to it.
…ing_focal_spot (cherry picked from commit a0f161c5d041fb1201fa976116cfccf2f2cb6600)
JevgenijaAksjonova
force-pushed
the
master
branch
from
3256e6c to
7655093
Compare
|
Hi! Are there more changes required? |
| __all__ = ('flying_focal_spot') | ||
|
|
||
|
|
||
| def flying_focal_spot(angle, apart, shifts): |
There was a problem hiding this comment.
Could you add a test for this in particular?
There was a problem hiding this comment.
What do you mean? Tests in geometry_tests use this function, since it is the only implemented option
There was a problem hiding this comment.
I agree, but they test the whole thing. It might be good to have a test for this in isolation to help track down errors.
There was a problem hiding this comment.
The other option, if this is not intended to be used externally, is to make this a private function in the cone beam geometry file.
|
Hi! Thank you for the feedback! |
|
Some of the seemingly outdated comments are still unresolved, could you check them and mark as resolved? |
|
Hi! Please check the latest updates. |
adler-j
left a comment
adler-j
left a comment
There was a problem hiding this comment.
Very minor comments, fix and go :)
| (the detector must be large enough, not to be influenced by shifts) | ||
| """ | ||
|
|
||
| # If no implementation is available, skip |
|
|
||
| # If no implementation is available, skip | ||
| if not odl.tomo.ASTRA_AVAILABLE: | ||
| pytest.skip(msg='ASTRA not available, skipping 2d test') |
There was a problem hiding this comment.
The fact that the test is skipped is documented automatically, I'd prefer information on why it was skipped. E.g. astra is needed to run it.
| assert np.max(proj[3, 15:]) > 5 | ||
|
|
||
|
|
||
| def test_source_detector_shifts_2d(): |
| phantom = odl.phantom.shepp_logan(space) | ||
|
|
||
| full_angle = np.pi | ||
| n_angles = 2 * 7 |
There was a problem hiding this comment.
it must be even
|
Hi, Jonas! The commit you checked was not the last one. I fixed some issues in ray_transform tests for 2d, now it works fine. im_combined =im_combined/ np.sum(im_combined) * np.sum(im) reduces the gap to 0.00079. Do you have any ideas, what is wrong? Could this be because the detector is not perpendicular to src_det_axis? |
|
How large is the relative error? 0.39 sounds like a lot (if the values should be ~1) so it might be worth investigating. Could you perhaps try to trigger it without the FFS somehow? |
JevgenijaAksjonova
force-pushed
the
master
branch
from
d91ba15 to
d53a35b
Compare
|
Hi! I have added detector shifts and changed the ray_transform tests so that source and detector positions coincide. Even after this a 2% error remained in cone beam, while using astra 1.8.3. Simply rescaling one of the images to match the scale of the other solves the issue. (This is not needed for the latest ASTRA, there it works fine). The detector shifts themself dont work as expected in 3D, this is related to issue #1567. |
|
Now I have added tests for the radon transform with detector shifts. It is ready for review |
|
Hi! Can I help with this somewhere? Also I have a question regarding the tam-danielsen window for ffs. To me it is not obvious how to mask out overlapping rays when the focal spot jumps between two different positions. Did you take care of this @JevgenijaAksjonova ? What way can this be solved? Maybe this was the source for the small error in 3d! |
adler-j
left a comment
adler-j
left a comment
There was a problem hiding this comment.
Some minor comments that should be doable without further review.
Sorry for the late review again.
|
|
||
| # If no implementation is available, skip | ||
| if not odl.tomo.ASTRA_AVAILABLE: | ||
| pytest.skip(msg='ASTRA not available, skipping 2d test') |
| det_shift_func=lambda angle: shift) | ||
|
|
||
| assert all_almost_equal(geom.angles, geom_shift.angles) | ||
| angles = geom.angles |
| geometries which mimic ffs by using initial angular offsets and | ||
| detector shifts | ||
| """ | ||
|
|
|
I was absent from GH for a while, sorry for the delay. Go ahead after the minor things @adler-j mentioned. There's nothing left to fix that can't be fixed afterwards :-) |
|
Okay, now I just did the nit edits myself. Merging after CI. |
|
Is the issue solved why @JevgenijaAksjonova got some errors in 3d? I think this has to do with the implementation of the tam-danielsen window. When we have projections at different source positions it is much more complicated to find and mask out the overlapping rays and this leads to errors in the backprojection. |
7 participants
Added possibility to arbitrary shift source positions. Since ASTRA allows to provide arbitrary source positions to projection and back-projection functions, it is possible now to do reconstruction with flying focal spot.
Implementation:
The sequence of unique shifts should be defined relative to the default source position and provided as an input to the geometry. Then this sequence is extended for all angles and returned when calling geometry.flying_focal_spot. However, one has to manually provide flying_focal_spot sequence to geometry.src_position() as an argument. This choice was made because FFS shifts correspond to discretization of angular partition and defining a FFS shift for any angle would introduce significant overhead in terms of work (moreover, in practice this is not needed).