Radon Transform of Discrete Data

[tknopp]

I am still exploring this package since I would like to reconstruct some real mirco-CT data.
In order to do so, I would like to forward-project the data that was reconstructed with the scanner software and compare it with the measured sinogram. However. It seams that radon is only defined for shape-phantoms like the sheep-logan phantom. Is there a way to pass an array of data and project that? In other words: Is there an interface like https://github.com/JuliaImages/ImageReconstruction.jl/blob/master/src/ImageReconstruction.jl#L15 that I am missing somehow?

[JeffFessler]

There currently is no efficient way to do that, but it is in the long-term plans and there is even a PR but it is WIP:
#14

There is an inefficient way though. In the same way that the Shepp-Logan phantom is just a (short) list of Ellipse objects, you could define a long list of Cuboid objects (one for each voxel), and then apply ImagePhantoms.radon to it. I suspect though that you will have millions of voxels and it might take a very long time. While you are waiting for it to run you could start Matlab and use the (mex compiled) forward projector I have over there and it might finish first, even counting the learning curve. Here's an example:
https://github.com/JeffFessler/mirt/blob/main/example/cone_beam_ct_example.m
Sadly the github version does not have the mex files so if you want to try it get the one from my web site:
https://web.eecs.umich.edu/~fessler/code/index.html

I will leave this issue open until we complete #14 which is what you really want.

[tknopp]

Thanks! Long-term plans is perfectly fine, I just want to stress test where this package is standing (geometry visualizations are absolutely wonderful! Want to integrate them in my lectures). A colleague of mine is currently using ASTRA and I want to see if it is possible to replace it with Sinograms.jl or not.

Regarding the inefficient way: This is certainly not an option because my sinogram is of size 2274x1013...

I have some other questions / discussion points (one regarding N(u)FFT based forward/backward projections and one regarding iterative reconstruction). Do you prefer discussions in Github/Issues or should I raise user questions in the discussion forum https://github.com/JuliaImageRecon/Sinograms.jl/discussions/landing?

[JeffFessler]

In the long run, having this be the Julia alternative to ASTRA is the goal. That is a high bar.
Meanwhile, I've initiated the discussion forum so anyone can post questions there.
If you need to make categories or such just go ahead (I am not to familiar with it yet).
Hopefully you have the ability to do so as an org member. If not let me know and I'll fix that.

[roflmaostc]

Hi,

as newcomer to the field, is that the algorithm you're talking about?
https://ieeexplore.ieee.org/abstract/document/1239600

Or related: https://ieeexplore.ieee.org/abstract/document/7866886

[JeffFessler]

is that the algorithm you're talking about?

yes, distance driven (or branchless version) is a wip here and there is PR #51 on it just waiting for me to have time to review it. (it will likely need quite a few changes.)

[roflmaostc]

Is there any news on it? @JeffFessler

I recently played with KernelAbstractions and had pretty good performance with generic code for CPU and GPU. Maybe I can polish the code and share it here. It's more like an educational radon transform based on interpolations, so quite naive. But performance is not that bad.

I wonder, if it would be worth to implement your draft with KernelAbstractions instead?

How do you do that in your lab? Do you have CUDA implementations? In CUDA C code?
I think, it's hard to implement an efficient radon transform in PyTorch, Jax, etc. because it's not easy (or impossible?) to write radon in a vectorized manner.

[JeffFessler]

I just merged an initial version of branchless distance driven (DD) projector / back-projector.
I won't tag this version because there are many updates I want to do with it, and the accuracy seems not quite as good as I would hope near the edges of the sinogram for some reason:
https://juliaimagerecon.github.io/Sinograms.jl/dev/generated/examples/08-bdd2d/
Our C code is CUDA based. It was painful. It would be great to have a KernelAbstractions version of any forward/back-projector pair! PRs welcome.

[roflmaostc]

I just published this: https://github.com/roflmaostc/RadonKA.jl

For me it really showed the power of Julia and KernelAbstractions. It's a very naive algorithm but I should compare it to your version here!

Especially I am interested in the exponential Radon transform (attenuation during travelling through the space.)

Maybe it is not too hard to change your code to KernelAbstractions? That would offer multithreading, CUDA, etc.
And your code is only 2D right? Generalization to 3D with KA would offer a good speedup instead of doing it slicewise, I believe.

[roflmaostc]

It's still not tagged, is it?

Would be interesting to test performance to RadonKA.jl

It also has a simple API to specify ray geometries.
I wanted to keep it simple and stupid. For my applications it seems to work fine.

[roflmaostc]

From the doc example:

Yours is fan beam whereas mine is parallel. But still a factor of 100 in difference.
Also your code seems to have some stability issues? Quite some allocations...

if !@isdefined(sinogramB)
    @time sinogramB = project_bdd(reverse(rot180(testimage'), dims=2), geo)
    sinogramB = sinogramB' # todo
end;
#   0.150224 seconds (220.91 k allocations: 217.402 MiB, 8.47% gc time)
# 
# and RadonKA.jl
@time sinogram_radonka =  RadonKA.radon(Unitful.ustrip(rot180(testimage)), π/2 .+ range(0, 2π, 90));
  0.001605 seconds (243 allocations: 241.156 KiB)

[JeffFessler]

I didn't tag yet because I would (eventually) like to update it a lot to improve allocations and type stability, but it will be a while before I can get to that.
If tagging the current version would be helpful let me know.