Tag Archives: Parallel projection

Radon and Hough Transform

Posted on by
3
Sinogram formed from Image:Sinogram Source - T...

Image via Wikipedia

This is to follow-up my earlier post.

In computed tomography of 2D parallel projection, its scanning process can be regarded as a real world example of 2D Radon transform.  Let us call density function (or beam attenuation coefficient or beam absorption coefficient) of body to be scanned f(x, y) .  Then scanned data for each projection angle \theta, after correcting exponential decay of beam due to attenuation, can be seen as Radon transform of f(x, y) :

\displaystyle Rf(\theta, X) = \int_{-\infty}^\infty f(X \cos \theta -Y \sin \theta, X \sin \theta + Y \cos \theta) dY.

Where (X, Y) are rotated coordinates:

\displaystyle x = X \cos \theta -Y \sin \theta,
\displaystyle y = X \sin \theta + Y \cos \theta.

Taking Fourier transform of f(x, y) we get:

\displaystyle F(\xi, \eta) = \int_{-\infty}^\infty \int_{-\infty}^\infty f(x, y) e^{-i(\xi x + \eta y)} dx dy.

Then F(\xi, \eta) in polar coordinates \xi = w \cos \theta, \eta = w \sin \theta would look like:

\displaystyle F(w \cos \theta, w \sin \theta) = \int_{-\infty}^\infty \int_{-\infty}^\infty f(x, y) e^{-i(\xi x + \eta y)} dx dy.

I want to change variables from (x, y) to (X, Y) in order to relate F to the Radon transform Rf.  Because of the change to polar coordinates, the exponent \xi x + \eta y becomes:

\displaystyle \xi x + \eta y = xw\cos \theta + yw \sin \theta = w(x \cos \theta + y \sin \theta) = wX.

And Jacobian is:

\displaystyle \dfrac{\partial x \partial y}{\partial X \partial Y} = \left | \begin{array}{ccc} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array} \right | = 1.

So F(w \cos \theta, w \sin \theta) is rewritten as:

\displaystyle F(w \cos \theta, w \sin \theta) = \int_{-\infty}^\infty [ \int_{-\infty}^\infty f(X \cos \theta - Y \sin \theta, X \sin \theta + Y \cos \theta) dY ] e^{-iwX} dX.

Gazing the inner integral with respect to Y for a while… carefully… This is just the projection of f in the direction of \theta, which is… Radon transform:

\displaystyle F(w \cos \theta, w \sin \theta) = \int_{-\infty}^\infty Rf(\theta, X) e^{-iwX} dX.

Okay, finally gazing this last integral for another while… patiently… It turns out that this is rightly Fourier transform of Rf with respect to X, doesn’t it?  My boss used to call this “Fourier slice theorem” or “central slice theorem“.  In CT, “filtered back projection” is a most important method for reconstructing f(x, y) from projection data Rf(\theta, X), and the method is based on the slice theorem.

But this theorem does not force us to go only backward (reconstruction) .  If we have an image f(x, y) at hand, then we can go forward to get lots of projection data using the same theorem.  Getting lots of projection data is same thing as plotting on (r, \theta) plane in Hough transform.  Steps to get projections using this idea are:

  1. Fourier transform the f(x, y) –> F,
  2. Resample F radially to convert to polar coordinates –> P,
  3. Inverse Fourier transform P –> projection data p.

Proper resampling should be difficult in practice.  It is like filtering in filtered back projection is difficult.  But I think it is still nice to be aware of the theorem, maybe just for fun.