computed tomography

Image via Wikipedia

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:

Fourier transform the $f(x, y)$ –> $F$ ,
Resample $F$ radially to convert to polar coordinates –> $P$ ,
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.

Got a book “Mathematical Methods in Image Reconstruction” by Frank Natterer and Frank Wübbeling. It brings about my old thought.

Say here is a picture of a box, or something rectangular. And let us assume that we roughly know position and attitude of the box before examining the picture. When we want to measure the position more accurately, sometimes we take the following procedure:

take orthogonal projection of (a certain interesting region of) the image to have 1 dimensional signal (integrate along parallel lines, in other words),
apply derivative filter (edge detector) to the 1D signal,
find peaks in the filtered signal.

Then we will know the position of 2 sides of the box.

But what if we do not know a priori the attitude of the box. We can not decide in which direction we should project the image.

One thing we can do for this case, if we still stick to the projection idea, is to project the image in many directions. And pick results that give the strongest edge filter response. Range of projection angle might be [0, 180) degrees if we know nothing about the pose of the box.

If we describe the procedure (up to projection) in mathematical way, we will reach 2D Radon transform.

In computer vision’s realm, there is a well-known technique for geometric feature extraction called Hough transform. In Hough transform for line detection, a most ineteresting step is to plot set of points in (r- $\theta$ ) plane for each feature point in the original image (x-y) plane. And I think this plotting step is essentially the same as the Radon transform.

I learned Hough transform when I was with computer vision company. But I did not know of Radon transform until I learned computer tomography things.

Some descriptions of Hough transform I met do not mention Radon transform, which is unfortunate because

0909

und sofort ins unendliche

Tag Archives: computed tomography

Radon and Hough Transform

Radon, Hough, projection to 1D