Reconstruction Algorithms and Scanning Geometries in Tomographic Imaging

the field of real numbers, and with $\rm R^2$ Euclidean plane, the elements of which are vectors: $\mathbf{r} \in \rm R^2$ means that $\mathbf{r}=(x,y)$ , where both x and y are real. For any two vectors a and b from $\rm R^2$ , their scalar product is defined by $\mathbf{a}\cdot \mathbf{b}={x_a}{x_b}+{y_a}{y_b}$ . With ${{{\mathbf{\uptau}}}_{\phi}}$ , we denote the unit vector ${{{ \mathbf{\uptau}}}_{\phi }} = (\cos \phi, \sin \phi)$ . The distribution of the specific parameter within the slice through the body is described by the object function $f : \rm D \to \rm R$ supported on the disk $D := \{(x,y); {x^2} + {y^2} \leq {R^2} \}$ . It is supposed that $f\in {\rm L_2}(\rm D)$ , which means that $\rm\int_D {{\it f^2}(\mathbf{r})\mathrm{ d}\mathbf{r}} <\infty,\ \mathbf{r}\in {R^2}$ . The Fourier transform of is defined by

$\hat{f}({ \mathbf{\uplambda}} )=\int_{\rm D} {f(\mathbf{r}){{\mathrm{ e}}^{{-{\rm i}(\mathbf{r}\cdot \boldsymbol{ \mathbf{\uplambda}} )}}}\mathrm{ d}\mathbf{r},} \quad \boldsymbol{ \mathbf{\uplambda}} \in {\rm R^2}.$

(6.1)

From the fact that the object function

is finitely supported, it follows that its Fourier transform is infinitely supported. Therefore, its inverse Fourier is expressed by

$f(\mathbf{r})=\frac{1}{{4{\pi^2}}}\int_{{{\rm R^2}}} {\hat{f}(\boldsymbol{ \uplambda} ){{\mathrm{ e}}^{{{\rm i}(\mathbf{r}\cdot \boldsymbol{ \mathbf{\uplambda}} )}}}\mathrm{ d}\boldsymbol{ \mathbf{\uplambda} }},$

(6.2)

where the integration is carried out over the whole plane ${\rm R^2}$ . The Radon transform of the function

can be defined by

$p(a,b)=\int_{{\rm D\cap \it L(a,b)}} {f(\mathbf{r}){\rm d}\mathbf{r}},$

(6.3)

where

is the line described by parameters

. The type of these parameters defines the type of the parameterization of the Radon transform. Commonly used is a $\phi t$ -parameterization:

$L(\phi, t):=\{(x,y):x\cos \phi +y\sin \phi =t\}.$

(6.4)

If parameter $\phi$ is fixed, then $L(\phi, t)$ generates a family of parallel lines. According to this, the parallel Radon projection in the direction $\phi$ is defined as

$p(\phi, t)=\int_{{\rm D\cap \it L(\phi, t)}} {f(\mathbf{r}){\rm d}\mathbf{r}} =\int_{{-\sqrt{{{R^2}-{t^2}}}}}^{{\sqrt{{{R^2}-{t^2}}}}} {f(t{{{\tau}}_{\phi }}+s\boldsymbol{ \tau}_{\phi}^{\bot })\mathrm{ d}s}, \quad -R\leq t\leq R$

(6.5)

(see Fig. 6.1, left).

Fig. 6.1

Parallel projection in direction $\phi$ (left); $\phi, t$ parameterization (right)

Depending on the situation, different notations, the sense of which will be always clear from the context, will be used, e.g., ${p_{\phi }}(t)$ , ${p_{\phi }}(f;t)$ , or $p(f;\phi, t)$ . If alternatively instead of t the angle parameter ψ is used (see Fig. 6.1, right), we speak of $\phi \psi$ parameterization. The line within this parameterization is expressed by

$x\cos \phi +y\sin \phi =R\cos \psi.$

(6.6)

Let ${x_{\alpha }}{y_{\alpha }}$ be coordinates in the coordinate system that is rotated relative to

at angle $\alpha$ :

$\begin{array}{lll} {x_{\alpha }}=x\cos \alpha +y\sin \alpha, \hfill \\{y_{\alpha }}=-x\sin \alpha +y\cos \alpha. \hfill\end{array}$

(6.7)

We speak of $\alpha \psi$ and $\alpha Y$ parameterizations of the Radon transform if

$L(\alpha, \psi ):=\{(x,y):{x_{\alpha }}\cos \psi +{y_{\alpha }}\sin \psi =R\cos \psi \}$

(6.8)

and

$L(\alpha, Y):=\{(x,y):Y{x_{\alpha }}+R{y_{\alpha }}=RY\},$

(6.9)

respectively. Families $L(\alpha, \psi )$ and $L(\alpha, Y)$ where only $\alpha$ is fixed generate fans of rays with an apex lying on the boundary of the disk $\mathrm{ D}$ . Therefore, both (6.8) and (6.9) are appropriate for the description of data collected over rays emitted by a point source moving along the boundary of the disk (Fig. 6.2).

Fig. 6.2

Fan-beam projection (left); setup for fan-beam geometry (right)

One of the important problems in practical tomography is the identification of optimal sampling conditions that must be imposed on the Radon transform. Because of limitations of the data acquisition system of the scanner, only discretized approximation of the Radon transform can be obtained. Therefore, on the one side, it is important to know how many data have to be collected in order to fully exhaust the potential of the scanner. On the other side, the radiation dose must be as small as possible, and therefore, the sparsest sampling scheme has to be identified among all of those which ensure the best approximation of the Radon transform. In connection to this problem, the Radon transform is supposed to be sampled on the lattice $\mathbf{L}$ , the elements of which are discrete points of the plane ${\rm R^2}$ coordinated by parameters of the Radon transform (e.g., parameters $\phi, t$ ). So-called admissible lattice can be generated by a $2\times 2$ matrix $\mathbf{W}$ such that $\mathbf{L}=\mathbf{W}{\rm Z^2}$ , i.e., for any $i,j\in \rm Z$ _,

${\mathbf{\upxi}} =\mathbf{W}\left( {\begin{array}{clclcllclcl}i \\j \\\end{array}} \right),\quad \boldsymbol{ \mathbf{\upxi}} \in \mathbf{L}.$

(6.10)

Lattice ${{\mathbf{L}}^{\bot }}$ is called dual or reciprocal to the lattice $\mathbf{L}$ if for any $\boldsymbol{ \mathbf{\upxi}} \in \mathbf{L}$ andany $\boldsymbol{ \mathbf{\upeta}} \in {{\mathbf{L}}^{\bot }}$ _, the scalar product $(\boldsymbol{ \mathbf{\upxi}}, {\mathbf{\upeta}} )=2\pi k,\ k\in \rm Z$ . The generator matrix of the reciprocal lattice is

${{\mathbf{W}}^{\bot }}=2\pi {{({{\mathbf{W}}^{-1 }})}^T},$

(6.11)

where the superscript

means transposition (see, e.g., [2]). Let the power spectrum of a function g be supported on the compact region

. The fundamental result of the sampling theory is as follows: if regions

and $K+{\mathbf{\upeta}}$ are disjoint for any $\boldsymbol{ \upeta} \in {{\mathbf{L}}^{\bot }}$ , then $g(\mathbf{r})=\sum\nolimits_{{\boldsymbol{ \mathbf{\upxi}} \in \mathbf{L}}} {g(\boldsymbol{ \mathbf{\upxi}} ){{\tilde{\chi}}_K}(\mathbf{r}-\boldsymbol{ \mathbf{\upxi}} )}$ , where ${{\tilde{\chi}}_K}$ is the inverse Fourier transform of the function ${\chi_K}$ , ${\chi_K}(\boldsymbol{ \mathbf{\uplambda}} )=1$ if $\boldsymbol{ \mathbf{\uplambda}} \in K$ , and otherwise ${\chi_K}(\boldsymbol{ \mathbf{\uplambda}} )=0$ (see [3]). In other words, the band-limited function can be exactly recovered from its discrete samples if the sampling conditions stated above are satisfied. From the definition of the reciprocal lattice, it follows that the denser ${{\mathbf{L}}^{\bot }}$ is, the sparser is $\mathbf{L}$ and vice versa. Therefore, if we want the sampling lattice $\mathbf{L}$ to be as sparse as possible, we have to find such a reciprocal lattice ${{\mathbf{L}}^{\bot }}$ that allows us to pack R₂ with mutually disjoint replicas of the

-region as tightly as possible. Certainly, the crucial role in doing this plays the shape of the

-region.

In some parts of the text, concepts from the theory of distributions (see, e.g., [4, 5]) are used. Within the framework of this theory, one defines the set $\mathrm{ B}$ of test functions. The distribution f over this set is to be understood as a map $f:\mathrm{ B}\to \rm R$ defined as follows: for $h\in \mathrm{ B}$ ,

$f[h]=\left\langle {f,h} \right\rangle =\int {f(x)\overline{h(x)}\mathrm{ d}x}.$

(6.12)

In the imaging science, such an approach is justified, e.g., when the output of the imaging system can be described in terms of the convolution with the impulse response $\hat{\varphi}$ . Let the function $\varphi :\mathrm{ R}\to \mathrm{ R}$ be smooth, symmetric, and compactly supported, and let its Fourier transform $\hat{\varphi}$ be essentially supported¹ on the interval

. Beyond this interval, $\hat{\varphi}(\omega )$ is required to decay faster than any degree of $|\omega {|^{-1 }}$ . We define the set $\mathrm{ B}$ as a family $\{{\varphi_t},t\in \rm R\}$ of translates of the impulse response e, i.e., ${\varphi_t}(s)=\varphi (s-t)$ . The set of Fourier transforms of test functions is $\bar{\mathrm{ B}} =\{{{\hat{\varphi}}_t},t\in \rm R\}$ , where

${{\hat{\varphi}}_t}(\omega )=\hat{\varphi}(\omega )\exp (-\mathrm{ i}\omega t).$

(6.13)

Then the Fourier transform of the distribution is defined as a map $\hat{f}:\bar{\mathrm{ B}} \to \rm R$ such that $\left\langle {\hat{f},{{\hat{h}}_t}} \right\rangle =2\pi \left\langle {f,{h_t}} \right\rangle$ for any $t\in \rm R$ .

6.3 Convolution-Based Algorithms

6.3.1 Parallel-Beam Geometry

The $\phi t$ parameterization is the most clear among all useful parameterizations of the Radon transform mentioned above. Being intensively investigated on the eve of practical tomography, the inversion of the Radon data collected within the parallel-beam geometry is fully understood today. Some important notions, such as, e.g., the resolution of the reconstruction, have been introduced and described for this kind of data geometry. Also the practical fan-beam inversion formulae are the result of the re-parameterization of the inversion formula obtained for the parallel-beam geometry.

Let ${p_{\phi }}(t)$ be the Radon projection of the object function $$ f $$

measured at angle $\phi$ . Due to (6.5), its 1D Fourier transform with respect to t is

${{\hat{p}}_{\phi }}(\omega )=\int_{-R}^R {\int_{{-\sqrt{{{R^2}-{t^2}}}}}^{{\sqrt{{{R^2}-{t^2}}}}} {f(t{{{\mathbf{\uptau}}}_{\phi }}+s{\mathbf{\uptau}}_{\phi}^{\bot }){{\mathrm{ e}}^{{-\mathrm{ i}\omega t}}}{\rm d}s{\rm d}t} }.$

(6.14)

Change of variables $\mathbf{r}=t{{{\mathbf{\uptau}}}_{\phi }}+s{\mathbf{\uptau}}_{\phi}^{\bot }$ , where $$ r=(x,y) $$

, in the double integral of (6.14) yields

${{\hat{p}}_{\phi }}(\omega )=\hat{f}(\omega \cos \phi, \omega \sin \phi ).$

(6.15)

The identity (6.15) is the statement of the so-called Fourier slice theorem. By taking the inverse Fourier transform, this identity can be rearranged to

$f(x,y)=\frac{1}{{8{\pi^2}}}\int_0^{{2\pi }} {\int_{{-\infty}}^{\infty } {\left| \omega \right|{{\hat{p}}_{\phi }}(\omega ){{{e}}^{{\mathrm{ i}\omega (\mathbf{r}\cdot {{\boldsymbol{ \uptau}}_{\phi }})}}}{\rm d}\omega {\rm d}\phi } }$

(6.16)

(see [6] for details). The formula (6.16) makes little sense if the product $|\omega |{{\hat{p}}_{\phi }}(\omega )$ is not absolutely integrable, what can be, e.g., in the case if ${{\hat{p}}_{\phi }}(\omega ){\sim} {\mkern6mu}{\omega^{{-2+\varepsilon }}},\varepsilon \geq 0$ . In order to stabilize it, let us multiply both sides of (6.15) with the term $|\omega |$ and handle the resulting identity in terms of distributions over the set B defined above. That is,

$\left\langle {|\omega |{{\hat{p}}_{\phi }}(\omega ),{{\hat{\varphi}}_t}(\omega )} \right\rangle =\left\langle {|\omega |\hat{f}(\omega \cos \phi, \omega \sin \phi ),{{\hat{\varphi}}_t}(\omega )} \right\rangle, \quad t\in \rm R$

(6.17)

or, written in explicit form,

$\int_{{-\infty}}^{\infty } {|\omega |{{\hat{p}}_{\phi }}(\omega )\hat{\varphi}(\omega ){{\mathrm{ e}}^{{\mathrm{ i}\omega t}}}{\rm d}\omega } =\int_{{-\infty}}^{\infty } {|\omega |\hat{f}(\omega \cos \phi, \omega \sin \phi )\hat{\varphi}(\omega ){{\mathrm{ e}}^{{\mathrm{ i}\omega t}}}\mathrm{ d}\omega }, \quad t\in \rm R.$

(6.18)

Setting $t=x\cos \phi +y\sin \phi$ , integrating over $\phi$ from $$ 0 $$

to $2\pi$ , rearranging the right-hand side appropriately, and changing the variables from polar to rectangular in the double integral on the right-hand side of (6.18) yield

$\frac{1}{2}\int_0^{{2\pi }} {\int_{{-\infty}}^{\infty } {|\omega |{{\hat{p}}_{\phi }}(\omega )\hat{\varphi}(\omega ){{\mathrm{ e}}^{{\mathrm{ i}\omega (\mathbf{r}\cdot {{\boldsymbol{ \uptau}}_{\phi }})}}}{\rm d}\omega \mathrm{ d}\phi } } =\int_{{{\rm R^2}}} {\hat{f}(\boldsymbol{ \mathbf{\uplambda}} )\hat{\varphi}(|\boldsymbol{ \mathbf{\uplambda}} |){{\mathrm{ e}}^{{\mathrm{ i}\left( {\mathbf{r}\cdot \boldsymbol{ \mathbf{\uplambda}}} \right)}}}{\rm d}\boldsymbol{ \mathbf{\uplambda} }}.$

(6.19)

The right integral in (6.19) is the inverse Fourier transform (up to the term $1/4{\pi^2}$ ) of the 2D convolution $$ f*E $$

, where

is a radial function such that $E(\mathbf{r})=\varphi (|\mathbf{r}|)$ . Hence, we obtain the inversion formula

$F(x,y)=\frac{1}{{8{\pi^2}}}\int_0^{{2\pi }} {\int_{{-\infty}}^{\infty } {|\omega |{{\hat{P}}_{\phi }}(\omega ){{\mathrm{ e}}^{{\mathrm{ i}\omega (\mathbf{r}\cdot {{\boldsymbol{ \uptau}}_{\phi }})}}}\mathrm{ d}\omega \mathrm{ d}\phi } },$

(6.20)

which has the same structure as (6.16), but

(6.21)

${P_{\phi }}(t)={p_{\phi }}*\varphi (t).$

(6.22)

Owing to the conditions imposed on $\varphi$ (see (6.13) and related discussion in Sect. 6.1), the formula (6.20) makes sense for any object function. Clearly, in terms of convolution, the formula (6.20) takes a form given by

$F(x,y)=\frac{1}{{4\pi }}\int_0^{{2\pi }} {\int_{-R}^R {P(\phi, t)H({t}^{\prime}-t)\mathrm{ d}t\mathrm{ d}\phi } },$

(6.23)

where ${t}^{\prime}$ is the parameter of the ray that goes through the point $$ (x,y) $$

and

$H(t)=\frac{1}{{2\pi }}\int_{-w}^w {|\omega |{{\mathrm{ e}}^{{\mathrm{ i}\omega t}}}\mathrm{ d}\omega }.$

(6.24)

The formula (6.23) is called filtered back projection (FBP) because the evaluation of $$ F $$

can be accomplished in two steps: the convolution of Radon projections with the kernel $$ H $$

(filtering) and evaluating $$ F(x,y) $$

with averages of ${{\{{P_{\phi }}*H(x\cos \phi +y\sin \phi )\}}_{{0\leq \phi < 2\pi }}}$ (back projection).

Note that (6.22), where on the right there is a convolution of parallel Radon projections with the impulse response of the imaging system, can be thought of as a model of measuring process. It is easy to check that

${P_{\phi }}(f;t)={p_{\phi }}(F;t),$

(6.25)

that is, the data obtained with the acquisition system of scanner are equivalent to Radon transform of the band-limited function F defined by (6.21). In [7], it was shown that the Fourier transform of the Radon transform of a band-limited function is essentially supported on the region, which has the shape of bow tie shown in Fig. 6.3, left. The rigorous estimate of this region is given in [8], p. 71:

$K=\{(k,\omega ):|\omega |< w,|k|<\max (|\omega |/\vartheta, (1/\vartheta -1)w)\},$

(6.26)

where $0<\vartheta <1$ is a parameter that is normally close to 1.

Fig. 6.3

Left:

-region of $p(F;\phi, t)$ ; w is the frequency bound and ${w_0}=(1/\vartheta -1)w$ . Middle: replicas of bow-tie K-region on the interlaced lattice $\mathbf{L}_i^{\bot }$ . Right: replicas of bow-tie $$ K $$

-region on the standard lattice $\mathbf{L}_s^{\bot }$

Following the sampling theory, let us define the lattice $$ L $$

$\mathbf{L}=\{({\phi_v},{t_{jv }}):{\phi_v}=2\pi v/N,{t_{jv }}=(j+vq/N)\Delta; \quad 0\leq q,v< N,\quad j\in \rm Z\},$

(6.27)

where

is supposed to be even, i.e., $$ N=2p $$

. As it follows from this definition, there are $$ N $$

directions ${\phi_v}$ ; for each direction, integrals over an equidistant set of parallel lines with spacing $\Delta$ are measured; the collection of parallel lines is shifted by an amount $$ vq/N $$

which varies with the angle ${\phi_v}$ . As mentioned above, $\mathbf{L}=\mathbf{W}{\rm Z^2}$ where $\mathbf{W}$ is a generator matrix and ${\rm Z^2}=\{(v,j),v,j\in \rm Z\}$ . From (6.27) it follows that

$\mathbf{W}=\left( {\begin{array}{clclcllclcl}{2\pi /N} & 0 \\{q\Delta /N} & \Delta \\\end{array}} \right).$

(6.28)

Then, the generator matrix for the reciprocal lattice is

${{\mathbf{W}}^{\bot }}=\left( {\begin{array}{clclcllclcl}N & {-q} \\0 & {2\pi /\Delta } \\\end{array}} \right)$

(6.29)

[see (6.11)]. Therefore, for the reciprocal lattice ${{\mathbf{L}}^{\bot }}$ _, we have the following representation:

${{\mathbf{L}}^{\bot }}=\left\{ {\left( {Nm-qk,\frac{{2\pi k}}{\Delta }} \right);\quad m,k\in \rm Z} \right\}.$

(6.30)

We will discuss two sampling lattices, namely, the standard lattice ${{\mathbf{L}}_s}$ ( $$ q=0 $$

) and the interlaced lattice ${{\mathbf{L}}_i}$ ( $$ q=p $$

). In the case of the standard sampling, the replicas of the $$ K $$

-region are mutually disjoint if

$p\geq w/\vartheta, \quad \Delta \leq \pi /w$

(6.31)

(see Fig. 6.3, right). Similarly, sampling conditions for the interlaced lattice are

$p\geq w(2-\vartheta )/\vartheta, \quad \Delta \leq 2\pi /w$

(6.32)

(see Fig. 6.3, middle). From (6.31) and (6.32), we see that the interlaced sampling lattice ${{\mathbf{L}}_i}$ is almost twice as sparser as the standard lattice ${{\mathbf{L}}_s}$ . If $\vartheta =1$ , then the interlaced sampling lattice can be obtained from the standard sampling lattice simply by dropping those grid points for which $$ v+j $$

is odd (Fig. 6.4). In other words, the same resolution as that of the data measured on the lattice ${{\mathbf{L}}_s}$ can be obtained with only one-half of the data by using the parallel-interlaced scanning geometry. In order to reconstruct data collected over the interlaced lattice, one could first interpolate the missing data to a denser standard lattice and then use (6.23) discretized on this denser lattice. The second approach would be to reconstruct directly from interlaced data (see [9] for error estimates of the corresponding reconstruction).

Fig. 6.4

The interlaced lattice (black circles) as a subset of the standard lattice (all intersection points). Both lattices have the same theoretical resolution

Practical application of the parallel-beam geometry can be found, e.g., in SPECT (single-photon emission computer tomography). The parallel-beam geometry of data collected with the SPECT device is a result of the collimation that is required for the correct association of the measured data with lines of integration. It has to be noted that because of the poorness of SPECT data and due to attenuation of the emitted photons in the object on their way to the detector, the FBP algorithm happens to be inappropriate. Different algorithms, mainly from the class of iterative methods, have been developed for the reconstruction of SPECT data (see review [10]).

In transmission tomography, the parallel-beam data can be collected with a pencil beam scanner of the first generation. The corresponding scanning process is a combination of N linear translations of a source-detector system across the object, each made along the line of the predefined orientation. The main disadvantage of this process is that data have to be literally built together in a “sample by sample” fashion, which results in a very long acquisition time.

6.3.2 Fan-Beam Geometry

Collecting fan-beam projections (see Fig. 3.2, left) is preferable in computed tomography because of a much shorter acquisition time compared to that required for the collection of parallel-beam projections.

Thorough derivation of corresponding inversion formulae can be found in [11] and [12]. Below, with minor modifications, we reproduce the derivation shown in [6].

The setup corresponding to a fan-beam geometry is as follows: the source moves around the object in the plane $$ xy $$

on the boundary of the disk $\mathrm{ D}$ . The center of the disk coincides with the origin of the coordinate system. If $\alpha$ is the parameter describing the position of the source, then at $\alpha =0$ the source is located on the $$ x $$

-axis. The coordinate system ${x_{\alpha }}{y_{\alpha }}$ defined by (6.7) rotates together with the source.

As it was mentioned above, Radon data collected with the point source can be efficiently described with parameters either $\alpha Y$ or $\alpha \psi$ (see (6.8) and (6.9) and Fig. 6.2, right). The fan beam, the rays of which are parameterized with ψ, is called equiangular. If the rays within the fan beam are parameterized with Y, the fan beam is said to be equidistant. In both cases, in order to obtain the inversion formulae, the formula (6.23) has to be re-parameterized properly.

In order to derive the formula for the inversion of the equidistant fan-beam data, the variables $(\phi, t)$ in the double integral of (6.23) have to be changed for the variables $(\alpha, Y)$ according to the transformation:

$\phi =\alpha +{\cot^{-1 }}\left( {\frac{Y}{R}} \right)$