CT Projection Simulator-2

I added a 2D Gaussian object to the projection simulator described in my previous post. I also cleaned up the code.

The Gaussian projections are computed from a formula derived in Bracewell’s book Two Dimensional Imaging[1] (see p. 532). This is a great book and I highly recommend it. To derive the Radon transform of a general Gaussian, Bracewell starts with a circularly symmetric Gaussian. Using his notation where ℛ is the Radon transform operator and R is the displacement of the line from the origin

(1) ℛ[e^− πr²] = e^− πR²

He writes e^− πr² = e^{− π(x² + y²)} and then applies a general coordinate transformation theorem with u = αx and v = βy getting

ℛ[e^{− π(α²x² + β²y²)}] = (exp[ − πR²α²β² ⁄ (β²cos²(θ) + α²sin²(θ))])/(√(β²cos²(θ) + α²sin²(θ))).

I think using scaling variables that specify the width of the peak is more intuitive so I substituted a = 1 ⁄ α and b = 1 ⁄ β resulting in

ℛ[e^{− π(x² ⁄ a² + y² ⁄ b²)}] = ab(e^{− πR² ⁄ s_m(θ)²})/(s_m(θ))

where s_m(θ)² = a²cos²(θ) + b²sin²(θ).

This formula is similar to the formula for the chord length of an ellipse from my previous post so I could re-use much of the code. The result is included in the code for this post. the code includes an example that creates the projections for two Gaussian blobs and reconstructs it. The results are shown in Fig. 1↓

Figure 1 Reconstructed image of two Gaussian blobs created with CTProjSim.

Bracewell’s derivation reminds me of an interesting subset of reconstruction algorithms based on the circular harmonic transform. Eq. 1↑ shows that the exponential is its own Radon transform. That is, it is an eigenfunction of the transform. You can quibble that the transform is really different from the original function--it has different physical units for example--but nevertheless it has the same functional form. Ein-Gal[3] developed a family of other “eigenfunctions” that are not circularly symmetric. Another member, for example is √(π)re^− πr²e^iθ. This is closely related to the Circular Harmonic Transform, which has been widely studied. It was used by Cormack[2] for example in his paper that got him part of the Nobel prize for CT. But as far as I know, it has not resulted in much practical application. The convolution-backprojection algorithm of Bracewell and Riddle[4] is almost universally used. I say almost because there has been a lot of recent interest in iterative algorithms to produce images with reduced noise.

Last edited Sep. 4, 2011

Linking is allowed but reposting or mirroring is expressly forbidden.

References

[1] Ronald Newbold Bracewell: Two Dimensional Imaging. Prentice Hall, 1994.

[2] A. M Cormack: “Representation of a Function by Its Line Integrals, with Some Radiological Applications”, Journal of Applied Physics, pp. 2722—2727, 1963.

[3] M. Ein-Gal: “The shadow transform: An approach to cross sectional imaging”, NASA STI/Recon Technical Report N, pp. 22706, 1974.

[4] R. N. Bracewell, A.C. Riddle: “”, The Astrophysical Journal, pp. 427—434, 1967. URL http://ecee.colorado.edu/~ecen4532/spring2010/labs/bracewell.pdf.