The Radon Transform

(1)

Department of Mathematics and Statistics, Villanova University, Villanova, PA, USA

Electronic supplementary material

The online version of this chapter (doi:10.1007/978-3-319-22665-1_2) contains supplementary material, which is available to authorized users.

2.1 Definition

For a given function f defined in the plane, which may represent, for instance, the attenuation-coefficient function in a cross section of a sample, the fundamental question of image reconstruction calls on us to consider the value of the integral of f along a typical line $\,\ell_{t,\,\theta }$ . For each pair of values of t and $\,\theta$ , we will integrate f along a different line. Thus, we really have a new function on our hands, where the inputs are the values of t and $\,\theta \,$ and the output is the value of the integral of f along the corresponding line $\,\ell_{t,\,\theta }$ . But even more is going on than that because we also wish to apply this process to a whole variety of functions f. So really we start by selecting a function f. Then, once f has been selected, we get a corresponding function of t and $\,\theta$ . Schematically,

$\displaystyle{\mathrm{input}\ f\,\mapsto \mathrm{output}\ \left \{(t,\,\theta )\mapsto \int \limits _{\ell_{t,\,\theta }}f\,ds\right \}.}$

This multi-step process is called the Radon transform, named for the Austrian mathematician Johann Karl August Radon (1887–1956) who studied its properties. For the input f, we denote by $\,\mathcal{R}(f)\,$ the corresponding function of t and $\,\theta \,$ shown in the schematic. That is, we make the following definition.

Definition 2.1.

For a given function f, whose domain is the plane, the Radon transform of f is defined, for each pair of real numbers $(t,\,\theta )$ , by

$\displaystyle\begin{array}{rcl} \mathcal{R}f(t,\,\theta )&:=& \int \limits _{\ell_{t,\,\theta }}f\,ds \\ & =& \int _{s=-\infty }^{\infty }f(t\cos (\theta ) - s\sin (\theta ),\ t\sin (\theta ) + s\cos (\theta ))\,ds.{}\end{array}$

(2.1)

A few immediate observations are that (i) both f and $\mathcal{R}f$ are functions; (ii) f is a function of the Cartesian coordinates x and y while $\mathcal{R}f$ is a function of the polar coordinates t and $\theta$ ; (iii) for each choice of f, t, and $\theta$ , $\mathcal{R}f(t,\,\theta )$ is a number (the value of a definite integral); (iv) in the integral on the right, the variable of integration is s, while the values of t and $\theta$ are preselected and so should be treated as “constants” when evaluating the integral.

We can visualize the Radon transform of a given function by treating $\,\theta \,$ and t as rectangular coordinates. We depict the values of the Radon transform according to their brightness on a continuum of grey values, with the value 0 representing the color black, 0. 5 representing a neutral grey, and the value 1 representing white. Such a graph is called a sinogram and essentially depicts all of the data generated by the X-ray emission/detection machine for the given slice of the sample. The choice of the term sinogram is no doubt suggested by the symmetry $\mathcal{R}f(-t,\,\theta +\pi ) = \mathcal{R}f(t,\,\theta )$ as well as by the appearance of the graphs for some simple examples that we will explore next. Sinograms can also be portrayed in color, using an appropriate segment of the rainbow in place of the grey scale.

2.2 Examples

Example 2.2.

As a first example, suppose our patient has a small circular tumor with radius 0. 05 centered at the point (0, 1), and suppose the attenuation-coefficient function f has the constant value 10 there. Now take an arbitrary value of $\,\theta$ . In this case, when $\,t =\sin (\theta )$ , the line $\,\ell_{t,\theta }\,$ will pass through (0, 1) and make a diameter of the circular tumor. Since the diameter of this disc is 2 ⋅ (0. 05) = 0. 1, we get $\mathcal{R}f(\sin (\theta ),\theta ) = 10 \cdot (0.1) = 1$ . Moreover, for any given $\,\theta$ , the value of $\,\mathcal{R}f(t,\,\theta )\,$ will be zero except on the narrow band $\,\sin (\theta ) - 0.05 < t <\sin (\theta ) + 0.05$ .

Thus, as $\,\theta \,$ varies from 0 to π, the graph will show a narrow, brighter grey ribbon of width 0. 1 centered around $\,t =\sin (\theta )$ . In other words, the graph of $\,\mathcal{R}f\,$ in the $\,(\theta,\,t)\,$ plane will resemble the graph of the sine function. Similarly, the graph in the $\,(\theta,\,t)\,$ plane of the Radon transform of a small, bright disc located at (1, 0) will resemble the graph of the cosine function. Figure 2.1 shows the sinograms for these two bright discs. Perhaps these examples helped to motivate the use of the term sinogram for the graph of a Radon transform.

Fig. 2.1

The sinograms for small, bright circular discs centered at (0, 1) (left) and at (1, 0) (right) resemble the graphs of the sine and cosine functions, respectively.

In the previous example, we considered an attenuation-coefficient function that had a constant (nonzero) value on a finite region of the plane and the value 0 outside of that region. To attach some terminology to functions of this type, suppose $\Omega$ is some finite region in the plane and take $f_{\Omega }$ to be the function that has the value 1 at each point contained in $\Omega$ and the value 0 at each point not in $\Omega$ . This function $f_{\Omega }$ is known as the characteristic function, or the indicator function, of the region $\Omega$ . Thus, in our previous example, we looked at 10 times the characteristic function of a small disc.

When the attenuation-coefficient function is a characteristic function $\,f_{\Omega }$ , its Radon transform is particularly easy to comprehend. Indeed, along any line $\,\ell_{t,\theta }$ , the value of $\,f_{\Omega }\,$ will be 0 except when the line is passing through the region $\,\Omega$ , where the value is 1. Thus, for each pair of values of t and $\theta$ , the value of $\mathcal{R}f_{\Omega }(t,\,\theta )$ is equal to the length of the intersection of the line $\ell_{t,\,\theta }$ with the region $\Omega$ .

In general, the object we are testing may comprise a collection of “blobs” of various materials, each with its specific attenuation coefficient. As we shall soon see, the Radon transform of the entire collection will be a composite of the transforms of the separate blobs. Consequently, understanding this basic sort of example will play a central role in assessing the accuracy of various approaches to image reconstruction.

Example 2.3.

For instance, let’s take the region $\Omega$ to be the closed circular disc of radius R > 0 centered at the origin. Then the characteristic function of $\Omega$ is given by

$\displaystyle{f_{\Omega }(x,\,y):= \left \{\begin{array}{cc} 1&\quad \mbox{ if $x^{2} + y^{2} \leq R^{\,2}$,} \\ 0& \quad \mbox{ otherwise.} \end{array} \right.}$

If we choose a value of t such that | t | > R, then, regardless of the value of $\theta$ , the line $\ell_{t,\,\theta }$ will not intersect the disc $\Omega$ . Thus, $f_{\Omega } = 0$ at every point on $\ell_{t,\,\theta }$ , and, hence, $\mathcal{R}f_{\Omega }(t,\,\theta ) = 0$ for such t. On the other hand, if | t | ≤ R, then $\ell_{t,\,\theta }$ intersects $\Omega$ along the segment corresponding to the parameter values $-\sqrt{R^{2 } - t^{2}} \leq s \leq \sqrt{R^{2 } - t^{2}}$ in the standard parameterization for the line. The length of this segment is $2\sqrt{R^{2 } - t^{2}}$ . The value of $f_{\Omega }$ is 1 at points in this interval and 0 otherwise. Therefore, for such t, the value of $\mathcal{R}f_{\Omega }(t,\,\theta )$ is the same as the length of the segment, namely, $2\sqrt{R^{2 } - t^{2}}$ .

In summary, we have shown that

$$$\displaystyle{\mathcal{R}f_{\Omega }(t,\,\theta ) = \left \{\begin{array}{cc} 2\sqrt{R^{\,2 } - t^{2}} & \quad \mbox{ if $\vert t\vert \leq R$,} \\ 0 &\quad \mbox{ if $\vert t\vert > R$.} \end{array} \right.}$$ ” src=”/wp-content/uploads/2016/10/A183987_2_En_2_Chapter_Equc.gif”></DIV></DIV></DIV></DIV></DIV><br /> <DIV id=FPar4 class=$

Example 2.4.

Continuing with the theme of characteristic functions, take $\mathcal{S}$ to be the region enclosed by the square whose edges lie along the vertical lines x = ±1 and the horizontal lines y = ±1. This square is centered at the origin and has side length 2. Let $f_{\mathcal{S}}$ be the function that takes the value 1 at each point of $\mathcal{S}$ and the value 0 at each point not in $\mathcal{S}$ . That is, let

$\displaystyle{f_{\mathcal{S}}(x,y) = \left \{\begin{array}{cc} 1&\quad \mbox{ if $\max \{\vert x\vert,\,\vert y\vert \}\leq 1$,}\\ 0 & \quad \mbox{ otherwise.} \end{array} \right.}$

As discussed above, for each t and each $\theta$ , the corresponding value $\mathcal{R}f_{\mathcal{S}}(t,\,\theta )$ , of the Radon transform of $f_{\mathcal{S}}$ , will equal the length of the intersection of the line $\ell_{t,\,\theta }$ and the square region $\mathcal{S}$ . There are two values of $\theta$ for which this intersection length is easy to see: $\theta = 0$ and $\theta =\pi /2$ . In these cases, the lines $\ell_{t,\,\theta }$ are vertical or horizontal, respectively. Thus, for $\theta = 0$ or $\theta =\pi /2$ , we get $\mathcal{R}f(t,\,\theta ) = 2$ when − 1 ≤ t ≤ 1 and $\mathcal{R}f(t,\,\theta ) = 0$ when | t | > 1. In general, the function $\mathcal{R}f_{\mathcal{S}}(t,\,\theta )$ will be piecewise linear in t for each fixed value of $\theta$ . Figure 2.2 shows these cross sections for several values of $\theta$ .

Fig. 2.2

Cross sections at $\theta = 0$ , $\theta =\pi /9$ , $\theta =\pi /4$ , and $\theta =\pi /3$ of the sinogram of the characteristic function of the basic square.

Figure 2.3 shows the full sinogram for the function $f_{\mathcal{S}}$ . Note the symmetry inherited from that of the square region $\mathcal{S}$ . This example will play an important role in the image reconstruction algorithms examined in Chapter 9

Fig. 2.3

The sinogram for the characteristic function of the square $\mathcal{S} =\{ (x,\,y)\,:\, \vert x\vert \leq 1,\,\mathrm{and}\,\vert y\vert \leq 1\}$ .

Example 2.5.

For an example that is not a characteristic function, let f be the function defined by

$$$\displaystyle{ f(x,\,y):= \left \{\begin{array}{cc} 1 -\sqrt{x^{2 } + y^{2}} & \quad \mbox{ if $x^{2} + y^{2} \leq 1$,} \\ 0 &\quad \mbox{ if $x^{2} + y^{2} > 1$.} \end{array} \right. }$$ ” src=”/wp-content/uploads/2016/10/A183987_2_En_2_Chapter_Equ2.gif”></DIV></DIV><br /> <DIV class=EquationNumber>(2.2)</DIV></DIV></DIV></DIV><br /> <DIV class=Para>The graph of <SPAN class=EmphasisTypeItalic>f</SPAN> is a cone, shown in Figure <SPAN class=InternalRef><A href=$ 2.4. We have already observed that, on the line $\ell_{t,\,\theta }$ , we have

$\displaystyle{x^{2} + y^{2} = (t\cos (\theta ) - s\sin (\theta ))^{2} + (t\sin (\theta ) + s\cos (\theta ))^{2} = t^{2} + s^{2}.}$

It follows that, on the line $\ell_{t,\,\theta }$ , the function f is given by

$$$\displaystyle\begin{array}{rcl} & & f(t\cos (\theta ) - s\sin (\theta ),\ t\sin (\theta ) + s\cos (\theta )) \\ & &:= \left \{\begin{array}{cc} 1 -\sqrt{t^{2 } + s^{2}} & \quad \mbox{ if $t^{2} + s^{2} \leq 1$,} \\ 0 &\quad \mbox{ if $t^{2} + s^{2} > 1$.} \end{array} \right.{}\end{array}$$ ” src=”/wp-content/uploads/2016/10/A183987_2_En_2_Chapter_Equ3.gif”></DIV></DIV><br /> <DIV class=EquationNumber>(2.3)</DIV></DIV>From this, we see that the value of <SPAN id=IEq89 class=InlineEquation><IMG alt=$ depends only on t and not on $\theta$ and that $\,\mathcal{R}f(t,\,\theta ) = 0\,$ whenever | t | > 1. For a fixed value of t such that | t | ≤ 1, the condition t ² + s ² ≤ 1 will be satisfied provided that s ² ≤ 1 − t ². Thus, for any value of $\,\theta \,$ and for t such that | t | ≤ 1, we have

$\displaystyle\begin{array}{rcl} & & f(t\cos (\theta ) - s\sin (\theta ),\ t\sin (\theta ) + s\cos (\theta )) \\ & &:= \left \{\begin{array}{cc} 1 -\sqrt{t^{2 } + s^{2}} & \quad \mbox{ if $-\sqrt{1 - t^{2}} \leq s \leq \sqrt{1 - t^{2}}$,} \\ 0 & \quad \mbox{ otherwise;} \end{array} \right.{}\end{array}$

(2.4)

whence

$\displaystyle{ \int _{\ell_{t,\,\theta }}f\,ds =\int _{ s=-\sqrt{1-t^{2}}}^{\sqrt{1-t^{2}} }\left (1 -\sqrt{t^{2 } + s^{2}}\right )\,ds\,. }$

(2.5)

This integral requires a trigonometric substitution for its evaluation. Sparing the details for now, we have

$\displaystyle{ \int _{s=-\sqrt{1-t^{2}}}^{\sqrt{1-t^{2}} }\left (1 -\sqrt{t^{2 } + s^{2}}\right )\,ds = \sqrt{1 - t^{2}} -\frac{1} {2}t^{2}\ln \left (\frac{1 + \sqrt{1 - t^{2}}} {1 -\sqrt{1 - t^{2}}}\right ). }$

(2.6)

Fig. 2.4

The figure shows the cone defined in (2.2) and the graph of its Radon transform for any fixed value of $\theta$ .

In conclusion, we have shown that the Radon transform of this function f is given by

$$$\displaystyle{ \mathcal{R}f(t,\,\theta ):= \left \{\begin{array}{cc} \sqrt{1 - t^{2}} -\frac{1} {2}t^{2}\ln \left (\frac{1+\sqrt{1-t^{2}}} {1-\sqrt{1-t^{2}}} \right )&\quad \mbox{ if $ - 1 \leq t \leq 1$,} \\ 0 & \quad \mbox{ if $\vert t\vert > 1$.} \end{array} \right. }$$ ” src=”/wp-content/uploads/2016/10/A183987_2_En_2_Chapter_Equ7.gif”></DIV></DIV><br /> <DIV class=EquationNumber>(2.7)</DIV></DIV>In this case, where <SPAN id=IEq93 class=InlineEquation><IMG alt=$ is independent of $\theta$ , the value of $\,\mathcal{R}f(t,\,\theta )\,$ corresponds to the area under an appropriate vertical cross section of the cone defined by z = f(x, y). Several of these cross sections are visible in Figure 2.4.

2.3 Some properties of $\mathcal{R}$

Suppose that two functions f and g are both defined in the plane. Then so is the function f + g. Since the integral of a sum of two functions is equal to the sum of the integrals of the functions separately, it follows that we get, for every choice of t and $\theta$ ,

$\displaystyle\begin{array}{rcl} \mathcal{R}(f + g)(t,\,\theta )& =& \int _{s=-\infty }^{\infty }(f + g)(t\cos (\theta ) - s\sin (\theta ),\ t\sin (\theta ) + s\cos (\theta ))\,ds {}\\ & =& \int _{s=-\infty }^{\infty }\{f(t\cos (\theta ) - s\sin (\theta ),\ t\sin (\theta ) + s\cos (\theta )) {}\\ & & \phantom{+ + +} + g(t\cos (\theta ) - s\sin (\theta ),\ t\sin (\theta ) + s\cos (\theta ))\}\,ds {}\\ & =& \int _{s=-\infty }^{\infty }f(t\cos (\theta ) - s\sin (\theta ),\ t\sin (\theta ) + s\cos (\theta ))\,ds {}\\ & & \phantom{++} +\int _{ s=-\infty }^{\infty }g(t\cos (\theta ) - s\sin (\theta ),\ t\sin (\theta ) + s\cos (\theta ))\,ds {}\\ & =& \mathcal{R}f(t,\,\theta ) + \mathcal{R}g(t,\,\theta ). {}\\ \end{array}$

In other words, $\mathcal{R}(f + g) = \mathcal{R}f + \mathcal{R}g$ as functions.

Similarly, when a function is multiplied by a constant and then integrated, the result is the same as if the function were integrated first and then that value multiplied by the constant; i.e., $\int \alpha f =\alpha \int f$ . In the context of the Radon transform, this means that $\mathcal{R}(\alpha f) =\alpha \mathcal{R}f$ .

We now have proven the following proposition.

Proposition 2.6.

For two functions f and g and any constants α and β,

$\displaystyle{ \mathcal{R}(\alpha f +\beta g) =\alpha \mathcal{R}f +\beta \mathcal{R}g. }$

(2.8)

In the language of linear algebra, we say that the Radon transform is a linear transformation; that is, the Radon transform $\,\mathcal{R}\,$ maps a linear combination of functions to the same linear combination of the Radon transforms of the functions separately. We also express this property by saying that “ $\,\mathcal{R}\,$ preserves linear combinations.”

Example 2.7.

Consider the function

$\displaystyle{f(x,\,y):= \left \{\begin{array}{cc} 0.5& \quad \mbox{ if $x^{2} + y^{2} \leq 0.25$,} \\ 1.0&\quad \mbox{ if $0.25 < x^{2} + y^{2} \leq 1.0$,} \\ 0 & \mbox{ otherwise.} \end{array} \right.}$

This a linear combination of the characteristic functions of two circular discs. Namely, $f = f_{\Omega _{1}} - (0.5)f_{\Omega _{2}}$ , where $\,\Omega _{1}\,$ and $\,\Omega _{2}\,$ are the discs of radii 1 and 0. 5, respectively, centered at the origin.

Using property (2.8), along with the computation in Example 2.3, it follows that

$$$\displaystyle\begin{array}{rcl} \mathcal{R}f(t,\,\theta )& =& \mathcal{R}(f_{\Omega _{1}})(t,\,\theta ) - (0.5)\mathcal{R}(f_{\Omega _{2}})(t,\,\theta ) {}\\ & =& \left \{\begin{array}{cc} 2\sqrt{1 - t^{2}} -\sqrt{(0.25) - t^{2}} & \quad \mbox{ if $\vert t\vert \leq 0.5$,} \\ 2\sqrt{1 - t^{2}} & \quad \mbox{ if $(0.5) < \vert t\vert \leq 1$,} \\ 0 & \quad \mbox{ if $\vert t\vert > 1$.} \end{array} \right.{}\\ \end{array}$$ ” src=”/wp-content/uploads/2016/10/A183987_2_En_2_Chapter_Equ10.gif”></DIV></DIV></DIV></DIV><br /> <DIV class=Para>Figure <SPAN class=InternalRef><A href=$ 2.5 shows the graph of this attenuation-coefficient function alongside a graph of the cross section of its Radon transform corresponding to any fixed value of $\,\theta \,$ and − 1 ≤ t ≤ 1.

Fig. 2.5

The figure shows the attenuation function defined in Example 2.7 alongside the graph of its Radon transform for any fixed value of $\theta$ .

What happens to the Radon transform if we modify a function either by shifting it or by re-scaling it? That is, suppose we know the Radon transform of a function f, and now look at the functions g(x, y) = f(x − a, y − b), where a and b are some real numbers, and h(x, y) = f(cx, cy), where c > 0 is a positive scaling factor.

In the first case, we obtain the graph of g by shifting the graph of f by a units in the x direction and b units in the y direction. It follows that if we take any line $\ell_{t,\,\theta }$ and shift it by just the right amount, we will get a line $\ell_{\widehat{t},\,\theta }$ with the property that $\mathcal{R}g(\widehat{t},\,\theta ) = \mathcal{R}f(t,\,\theta )$ . What is the correct shift in the value of t? Well, when t = 0, then $\ell_{0,\,\theta }$ passes through the origin, while the parallel line $\ell_{a\cos (\theta )+b\sin (\theta ),\,\theta }$ passes through (a, b), as we saw in Exercise 12 of Chapter 1 So, the relationship between $\ell_{0,\,\theta }$ and the graph of f is the same as that between $\ell_{a\cos (\theta )+b\sin (\theta ),\,\theta }$ and the graph of g. In other words, the correct shift is $\widehat{t} = t + a\cos (\theta ) + b\sin (\theta )$ . Figure 2.6 illustrates this correspondence.

Fig. 2.6

The shifting property of $\mathcal{R}$ : The lines $\ell_{t,\theta }$ and $\ell_{t+a\cos (\theta )+b\sin (\theta ),\theta }$ intersect congruent regions, centered at (0, 0) and (a, b), respectively, along segments of the same length.

In the case of the function h(x, y) = f(cx, cy), we can think of the domain of h as a (1∕c)-times scale model of the corresponding domain for f. Thus, to compute the Radon transform of h for a given choice of t and $\theta$ , we first have to scale the value of t by the factor c, in order to locate a parallel line that intersects a similar part of the domain of f. Then multiply the corresponding value of the Radon transform of f by 1∕c to get back to the scale of h. That is, for given values of t and $\theta$ , we get $\mathcal{R}h(t,\,\theta ) = (1/c) \cdot \mathcal{R}f(ct,\,\theta )$ . This relationship is shown in Figure 2.7, in which f and h are taken to be the characteristic functions of two square regions.