If you have trouble viewing this, try the pdf of this post. You can download the code used to produce the figures in this post.

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
He writes *e*^{ − πr2} = *e*^{ − π(x2 + y2)} and then applies a general coordinate transformation theorem with *u* = *α**x* and *v* = *β**y* getting *a* = 1 ⁄ *α* and *b* = 1 ⁄ *β* resulting in

ℛ[*e*^{ − π(α2x2 + β2y2)}] = (exp[ − *π**R*^{2}*α*^{2}*β*^{2} ⁄ (*β*^{2}cos^{2}(*θ*) + *α*^{2}sin^{2}(*θ*))])/(√(*β*^{2}cos^{2}(*θ*) + *α*^{2}sin^{2}(*θ*))).

I think using scaling variables that specify the width of the peak is more intuitive so I substituted
ℛ[*e*^{ − π(x2 ⁄ a2 + y2 ⁄ b2)}] = *ab*(*e*^{ − πR2 ⁄ sm(θ)2})/(*s*_{m}(*θ*))

where *s*_{m}(*θ*)^{2} = *a*^{2}cos^{2}(*θ*) + *b*^{2}sin^{2}(*θ*).

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↓

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*^{ − πr2}*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

© 2011 by Aprend Technology and Robert E. Alvarez

Linking is allowed but reposting or mirroring is expressly forbidden.

[1] : *Two Dimensional Imaging*. Prentice Hall, 1994.

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

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

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