JuliaImageRecon / Sinograms.jl

Julia library for working with sinograms / tomography / Radon transform
MIT License
15 stars 6 forks source link

Parker Weights: Fan Beam - Short Angle #39

Closed tknopp closed 1 year ago

tknopp commented 1 year ago

My excitement for this package is rapidly growing in the last days. Today I had a breakthrough in reconstructing experimental flat-panel fan-beam data from this scanner https://www.scanco.ch/xtremectii.html after understanding/tweaking the geometry. The Sinograms.jl reconstruction now matches the reconstruction from the vendor.

On the way I required implementing the parker weights for short angle fan beam data. I did this outside Sinograms.jl in this way:

plan = plan_fbp(sg, ig, window=Window(Hamming(), filter))
plan.parker_weight .= 1.0
parkerWeights = calcParkerWeights(R)
fbp_image, sino_filt = fbp(plan, R.*parkerWeights)

with

function calcParkerWeights(R)
  parkerWeights = ones(size(R))
  λ = range(0,R.angleRange, size(R,2))
  γ = range(-R.fanAngle/2,R.fanAngle/2, size(R,1))
  γThresh = R.fanAngle/2 
  for y=1:size(R,2)
      for x=1:size(R,1)
        if 0 <= λ[y] < 2*(γThresh + γ[x])
            parkerWeights[x,y] = sin(π/4* λ[y]/(γThresh + γ[x]) )^2
        end
        if λ[end] - 2*γThresh + 2*γ[x] < λ[y] <= λ[end] && γThresh + γ[x] != 0.0
            parkerWeights[x,y] = sin(π/4* (λ[end]-λ[y])/(γThresh - γ[x]) )^2
        end       
      end
  end
  reverse(parkerWeights,dims=1)
end

Here, R is an ImageMeta object and has some additional metadata. When integration the code into Sinograms.jl we would of course take the necessary information from the Sinogram geometry object. The mathematics of the code is based on https://opus4.kobv.de/opus4-fau/files/801/dissertation.pdf (page 26, formula 3.48) with a minor correction (minus sign).

Now the question: How can we integrate this in Sinograms.jl? I did not study the source code and thus am not totally sure. What irritates me is that plan.parker_weight is 1D while the parker weights I implemented from the reference are 2D. Any hint?

JeffFessler commented 1 year ago

The 1D Parker weights was a bug. (I just had not gotten around to implementing it yet. All it took was someone being interested 😄) PR #40 addresses this and includes a translation of my Matlab code into Julia. I didn't use your code because pretty much everything here is a translation of Matlab to Julia so I was confident it would work and it was already using my variable names. Now we can compare yours vs mine and see if there are differences that matter. After the CI passes, would you mind trying it out?

JeffFessler commented 1 year ago

The tests passed, with fewer iterations than I expected... I can wait for your code review, or just merge if you prefer and then wait to tag it until you test drive it and give an OK. Thanks for the nudge!

tknopp commented 1 year ago

My local tests on the experimental data worked fine as well. And yes the intention was not to preserve anything from my code. It was just meant to showcase what I have done.

JeffFessler commented 1 year ago

Great! I will tag after some further tweaks in #41 pass.