NuFFT based Forward/Backward Projections

[tknopp]

I was wondering if N(u)FFT based Forward/Backward Projections are in the scope of this package (e.g. https://ieeexplore.ieee.org/document/1398550). I don't have too much practical experience with experimental CT data so that I can't judge why they are not used so much but my understanding is that they are faster than FBP (N^2 log N^2 instead of N^3 for NxN images). So is there some limitation (e.g. flexibility in the geometry?) which makes the N(u)FFT approach not fit well into Sinograms.jl?

[JeffFessler]

They could fit in, certainly. Indeed I have 2D versions in Matlab MIRT that could be translated. But AFAIK it only works natively for 2D parallel beam and some interpolation is needed to do 2D fan beam and I think most (all?) "real" CT work is 3D these days so it would be a low priority for me I think. Happy to hear more...

[tknopp]

Thanks! Good to know that in principle it could fit here. The linked paper does apply it to fan beam data so I wonder if interpolation is really a concern here.

In any case it certainly makes sense to first have some space-based forward projector (working on array data) so that one can implement iterative solvers based on that.

[roflmaostc]

Hi,

in best case a forward projector would be also ADable. Then simple inclusion in a ML library would be possible.

Is there anything planned in that direction?

Best,

Felix

[JeffFessler]

Forward and backprojection are linear so using raw AD seems unlikely to be efficient. We simply use one line of ZygoteRules.@adjoint in the SPECT case here:
https://juliaimagerecon.github.io/SPECTrecon.jl/stable/generated/examples/6-dl/#Custom-backpropagation
and anyone can do the same for the CT case.

[roflmaostc]

Ah, makes sense.
Is the (unfiltered?) backprojection always the adjoint? Also in cone beam geometries?

Can you recommend a reference where that's explained in more detail?

[JeffFessler]

As long as one implements the back-projector to be the adjoint of the forward projector then @adjoint will be appropriate by definition. It takes effort to verify correctness. Any code that I share will have such checks, e.g.:
https://juliaimagerecon.github.io/SPECTrecon.jl/stable/generated/examples/2-rotate/#Adjoint
https://github.com/JuliaImageRecon/SPECTrecon.jl/blob/main/test/adjoint-project.jl

I'm not sure of a reference that discusses the equality <A*x,y> = <x,A'y> that is the foundation for some of these checks...
But this reference is useful:
https://doi.org/10.1088/0031-9155/49/11/024