Singular-value decomposition of a tomosynthesis system

Anna Burvall; Harrison H. Barrett; Kyle J. Myers; Christopher Dainty

doi:10.1364/OE.18.020699

1. Introduction

The tomosynthesis method [1, 2] is currently considered to be of great potential, especially in detection of breast cancer [3]. The idea is to replace mammography, which gives two-dimensional (2D) projection, with tomosynthesis which gives three-dimensional (3D) information e.g. as slice images. The ultimate purpose is to improve the detection of breast cancer, at a relatively small increase in dose and cost. Prototype tomosynthesis systems have already been developed and pilot clinical trials performed [4–7]. These trials have so far not demonstrated distinct advantages in detection compared to mammography, but indicate that tomosynthesis images have better visualization of structures without masking by overlapping tissue. Improved detection has been clearly demonstrated in studies on mastectomy specimens [8].

Tomosynthesis can be said to be half-way between a single projection as used in mammography, and the full-angle scan performed in computed tomography (CT). The first produces a 2D representation and the other a 3D representation. Tomosynthesis, where a number of projections are taken at a limited angle, provides partial 3D information — the resolution in one direction must be limited. The hope is that tomosynthesis will provide the advantage of 3D information over conventional mammography, without the increased dose and cost of a full CT.

The 3D object, for example, the breast tissue, must be reconstructed from the 2D projection images. Several techniques are used [9,10], mainly back-projection methods, iterative methods, or statistical methods. The reconstruction task is more difficult than in tomography, since the information is incomplete and a true inverse transform cannot be found. Especially in analytical methods such as back-projection, this has a tendency to produce artifacts. We propose to use singular-value decomposition (SVD) [11–13] for the reconstruction, in a manner very similar to the technique we developed for through-focus imaging systems [14]. The result produced is the measurement component of the object, which could be regarded as its pseudoinverse [11] reconstruction and is as close to the inverse as we can get. The measurement component is the part of the object that is actually transmitted into the image. An advantage of the method is its mathematical stringency: what does not show up in the measurement component simply cannot be reconstructed, and the measurement component is unique for each set of measurements. A disadvantage is that the measurement component is different from a conventional reconstruction since it contains negative values. The measurement component could be used for further studies on image quality and detectability, using a channelized Hotelling observer [11, 15, 16].

At present, we are analyzing a parallel-path setup where the x-ray source moves in a plane parallel to the detector plane. Other geometries [2] could be analyzed but it would require some work to adjust the analysis. We consider the geometrical approximation, disregarding any effects of diffraction. The object is seen as purely absorbing, with no scattering or phase contribution. Scattering cannot be included, since it breaks the lateral shift-invariance. As we analyze the ideal case, we do not yet consider noise.

2. SVD analysis

We consider a geometry where all source positions are in one transverse plane, and images captured in a second plane parallel to the first and placed behind the object. This geometry is illustrated in Fig. 1. In the figure the object is shown as limited in both z, x and y direction. We present this case of limited object and image support as the most general, but will mainly consider the simplified case of infinite transverse extent. The support in the z direction is per definition limited to at most the distance between the source and image planes. The variables used are r = (x,y) for transverse coordinates in object space and z for longitudinal coordinate in object space. The image plane is placed at the origin (z = 0), and has transverse coordinates r _d = (x_d,y_d). The source plane is placed at z = −s and position number m given by r _m = (χ_m, ψ_m). The source is assumed to be a point. The object f(r,z) is identical to the absorption coefficient, the image is g_m(r _d) for the image taken using source position r _m, and g(r _d) is the vector of all images.

Fig. 1 Geometry of the considered tomographic system. The different source positions are χ_m, and the screen is placed along the z axis.

Download Full Size | PDF

A relation between object and image can be found from the steady-state solution of the Boltzmann transport equation in absorbing media, assuming a point source [11]:

I_{m} (r_{d}) = \frac{A cos θ}{{| r_{1} - r_{0} |}^{2}} exp [- \int_{0}^{| r_{1} - r_{0} |} d l' f (r_{1} - \frac{r_{1} - r_{0}}{| r_{1} - r_{0} |} l')],

where I_m(r _d) is the irradiance, A is the strength of the point source in emitted flux per unit volume per unit solid angle, r ₁ = (r _d, 0) = (x_d, y_d, 0) is the position of the image point in three dimensions and r ₀ = (r _m, −s) = (χ_m, ψ_m, −s) the position of the source. The angle θ is the angle of incidence of the ray from r _m to r _d onto the detector plane. A change of variables to our chosen coordinate system, namely z = −l′s/|r ₁ − r ₀|, leads to

log (\frac{I_{m} (r_{d}) {| r_{1} - r_{0} |}^{2}}{A cos θ}) = \frac{| r_{1} - r_{0} |}{s} \int_{0}^{- s} d z f (r_{0} + (r_{0} - r_{1}) \frac{z}{s}) .

Using cosθ = s/|r ₁ − r ₀| and |r ₁ − r ₀| = (|r _d − r _m|² + s ²)^1/2, and changing notation from the three-dimensional r ₁ and r ₀ to the two-dimensional r _d and r _m, leads to

log [\frac{1}{A s} {({| r_{d} - r_{m} |}^{2} + s^{2})}^{3 / 2} I_{m} (r_{d})] = \int_{0}^{- s} d z f (r_{d} + (r_{d} - r_{m}) \frac{z}{s}, z)

We introduce a relation between the irradiance I_m(r _d) registered in the image plane and the image function g_m(r _d):

g_{m} (r_{d}) = \frac{s}{{(s^{2} + {| r_{d} - r_{m} |}^{2})}^{1 / 2}} log [\frac{1}{A s} {(s^{2} + {| r_{d} - r_{m} |}^{2})}^{3 / 2} I_{m} (r_{d})] .

This definition of g_m(r _d) is designed to yield a linear relation between object and image, and inserting Eq. (2) into Eq. (1) leads to the propagation integral from object to image

g_{m} (r_{d}) = {[H f]}_{m} (r_{d}) = \int_{z_{1}}^{z_{2}} d z \iint_{- \infty}^{\infty} d^{2} r b_{f} (r, z) f (r, z) h_{m} (r, r_{d}, z),

and its adjoint (or backprojection) operator

[H^{†} g] (r, z) = \sum_{m = 1}^{M} \iint_{- \infty}^{\infty} d^{2} r_{d} b_{g} (r_{d}) g_{m} (r_{d}) h_{m}^{*} (r, r_{d}, z),

where z ₁ and z ₂ limit the object extension in the z direction, −s < z ₁ < z ₂ < 0, M is the number of source positions, b_f(r,z) and b_g(r _d) define the object and image support respectively, the asterisk denotes complex conjugate, and the impulse-response function h is the two-dimensional Dirac delta function

h_{m} (r, r_{d}, z) = δ (r - r_{d} - (r_{d} - r_{m}) \frac{z}{s}) .

The object and image support contained in b_f(r,z) and b_g(r _d) can in general be adjusted according to specific situations, but the support in the axial direction can never extend beyond the source or image plane as indicated by the limits z ₁ and z ₂. Eqs. (3)–(5) give a complete description of the propagation, and consequently also of the reconstruction. Now the SVD method for multiple images, earlier presented for telecentric optical imaging systems [14], can be applied. However, this method requires shift invariance in the transverse direction, so we apply a change of variables [17]:

\tilde{r} = (\tilde{x}, \tilde{y}) = (\frac{s}{s + z} x, \frac{s}{s + z} y), \tilde{z} = \frac{s}{s + z} z

with the inverse relations

r = (x, y) = (\frac{s}{s - \tilde{z}} \tilde{x}, \frac{s}{s - \tilde{z}} \tilde{y}), z = \frac{s}{s - \tilde{z}} \tilde{z} .

Intuitively, the transformation in Eq. (6) means the source plane is moved to minus infinity and the rays from one source point become parallel as illustrated in Fig. 2. The position r _m of the source determines the inclination of the rays. The object is distorted by this transformation, but the image remains the same. In these coordinates, the impulse-response function becomes

{\tilde{h}}_{m} (\tilde{r} - r_{d}, \tilde{z}) = {(\frac{s - \tilde{z}}{s})}^{2} δ (\tilde{r} - r_{d} + \frac{\tilde{z}}{s} r_{m}),

Fig. 2 The change of variables moves the source plane to minus infinity. Rays originating from two source positions, one on-axis and one off-axis, are shown. The object in this figure is distorted compared to the object in Fig. 1, but the image is not.

Download Full Size | PDF

where h̃_m(r̃ − r _d,z̃) = h_m(r,r _d,z). Since the Jacobian of the transformation is

| \frac{\partial (x, y, z)}{\partial (\tilde{x}, \tilde{y}, \tilde{z})} | = \frac{s^{4}}{{(s - \tilde{z})}^{4}}

the integral relations now become

g_{m} (r_{d}) = {[H \tilde{f}]}_{m} (r_{d}) = \int_{{\tilde{z}}_{1}}^{{\tilde{z}}_{2}} d \tilde{z} \iint_{- \infty}^{\infty} d^{2} \tilde{r} \frac{s^{4}}{{(s - \tilde{z})}^{4}} {\tilde{b}}_{f} (\tilde{r}, \tilde{z}) \tilde{f} (\tilde{r}, \tilde{z}) {\tilde{h}}_{m} (\tilde{r} - r_{d}, \tilde{z}),

and

[H^{†} g] (\tilde{r}, \tilde{z}) = \sum_{m = 1}^{M} \iint_{- \infty}^{\infty} d^{2} r_{d} b_{g} (r_{d}) g_{m} (r_{d}) {\tilde{h}}_{m}^{*} (\tilde{r} - r_{d}, \tilde{z}),

where f̃(r̃,z̃) = f(r,z) is the distorted object function and b̃_f(r̃,z̃) = b_f(r,z) the distorted object-support function. The integration limits are z̃ ₁ = sz ₁/(s + z ₁) and z̃ ₂ = sz ₂/(s + z ₂) and fulfill −∞ < z̃ ₁ < z̃ ₂ < 0. We note that due to the change of variable, the distorted object- and image-space functions are defined on the Hilbert spaces with weighted inner products [18]

{〈 {\tilde{f}}_{1} (\tilde{r}, \tilde{z}), {\tilde{f}}_{2} (\tilde{r}, \tilde{z}) 〉}_{obj} = \iint_{- \infty}^{\infty} d^{2} \tilde{r} \int_{{\tilde{z}}_{1}}^{{\tilde{z}}_{2}} d \tilde{z} {\tilde{b}}_{f} (\tilde{r}, \tilde{z}) {\tilde{f}}_{1} (\tilde{r}, \tilde{z}) {\tilde{f}}_{2}^{*} (\tilde{r}, \tilde{z}) \frac{s^{4}}{{(s - \tilde{z})}^{4}}

and

{〈 g_{1} (r_{d}), g_{2} (r_{d}) 〉}_{im} = \sum_{m = 1}^{M} \iint_{- \infty}^{\infty} d^{2} r_{d} b_{g} (r_{d}) g_{1 m} (r_{d}) g_{2 m}^{*} (r_{d}) .

The weighted inner product of Eq. (11) is required in order to fulfill the adjoint relation

{〈 H \tilde{f}, g 〉}_{im} = {〈 \tilde{f}, H^{†} g 〉}_{obj} .

Combining Eqs. (9) and (10) finally yields the image-space Hermitian operator as

{[H, H^{†} g]}_{m} (r_{d}) = \sum_{m' = 1}^{M} \iint_{- \infty}^{\infty} d^{2} r_{d}^{'} b_{g} (r_{d}^{'}) g_{m'} (r_{d}^{'}) k_{m m'} (r_{d}, r_{d}^{'})

where

k_{m m'} (r_{d}, r_{d}^{'}) = \int_{{\tilde{z}}_{1}}^{{\tilde{z}}_{2}} d \tilde{z} \iint_{- \infty}^{\infty} d^{2} \tilde{r} {\tilde{b}}_{f} (\tilde{r}, \tilde{z}) \frac{s^{4}}{{(s - \tilde{z})}^{4}} {\tilde{h}}_{m} (\tilde{r} - r_{d}, \tilde{z}) {\tilde{h}}_{m'}^{*} (\tilde{r} - r_{d}^{'}, \tilde{z}) .

Similarly, the object-space Hermitian operator will be

[H^{†} H \tilde{f}] (\tilde{r}, \tilde{z}) = \int_{{\tilde{z}}_{1}}^{{\tilde{z}}_{2}} d \tilde{z}' \iint_{- \infty}^{\infty} d^{2} \tilde{r}' {\tilde{b}}_{f} (\tilde{z}', \tilde{r}') \tilde{f} (\tilde{z}', \tilde{r}') p (\tilde{r}, \tilde{r}', \tilde{z}, \tilde{z}')

where

p (\tilde{r}, \tilde{r}', \tilde{z}, \tilde{z}') = \frac{s^{4}}{{(s - \tilde{z})}^{4}} \sum_{m = 1}^{M} \iint_{- \infty}^{\infty} d^{2} r_{d} b_{g} (r_{d}) {\tilde{h}}_{m}^{*} (\tilde{r} - r_{d}, \tilde{z}) {\tilde{h}}_{m} (\tilde{r}' - r_{d}, \tilde{z}') .

2.1. Infinite object and image

If both object and image are considered infinite in the transverse direction, and the object confined in the longitudinal direction by z ₁ < z < z ₂, both support functions b_f(r, z) and b_g(r _d) are identically equal to one. While this assumption is obviously not true, its effects on the resulting analysis are surprisingly small. Its main contribution is edge effects. Theoretically, the object and images are assumed to be infinite, but for the evaluation, finite object and images must be used. This leads to strange effects close to the edges, so the rim of the resulting images includes non-reliable results and should be discarded. The assumption simplifies the analysis considerably, since the system is now transversally shift-invariant and the kernels of the Hermitian operators may be written k _mm′(r _d − r ′_d) and p(r̃− r̃′,z̃,z̃′). Insertion of Eq. (8) into Eq. (13) yields, after performing the integral in r̃,

k_{m m'} (r_{d} - {r'}_{d}) = \int_{{\tilde{z}}_{1}}^{{\tilde{z}}_{2}} d \tilde{z} δ (r_{d} - {r'}_{d} - \frac{\tilde{z}}{s} (r_{m} - {r'}_{m})) .

This integral could easily be evaluated, but the subsequent analysis gets easier if it is left as it is. Assuming singular functions of the form

v_{ρ, j} (r_{d}) = V_{j} (ρ) exp (2 π i ρ \cdot r_{d})

yields, after insertion into Eq. (12), the eigenequation [14]

K (ρ) V_{j} (ρ) = μ ρ,_{j} V_{j} (ρ)

where $K$ (ρ) is an M × M matrix with elements

k_{m m'} = \iint_{- \infty}^{\infty} d^{2} r_{d} k_{m m'} (r_{d}) exp (- 2 π i ρ \cdot r_{d}) .

Inserting Eq. (14) into Eq. (17) and performing the integral in r _d yields

k_{m m'} = \int_{{\tilde{z}}_{1}}^{{\tilde{z}}_{2}} d \tilde{z} exp (- 2 π i ρ \cdot (r_{m} - r_{m'}) \frac{\tilde{z}}{s}),

which is integrated to yield

K_{m m'} = {\tilde{z}}_{2} - {\tilde{z}}_{1}

for m = m′ or ρ = (0,0), and

K_{m m'} = \frac{i s}{2 π ρ \cdot (r_{m} - r m')} {exp [- 2 π i ρ \cdot (r_{m} - r m') \frac{{\tilde{z}}_{2}}{s}] - exp [- 2 π i ρ \cdot (r_{m} - r m') \frac{{\tilde{z}}_{1}}{s}]}

for m ≠ m′. Thus the elements of $K$ have been found, and Eq. (16) can be solved numerically to yield the vectors V _j(ρ). The image-space singular functions have been found, in exactly the same manner as Burvall et al. [14]. The object-space eigenfunctions can then be found by projection of the image-space eigenfunctions using the adjoint operator in Eq. (10).

The solution can also be obtained in object space. The results will be identical to those obtained from the image-space analysis, but we still outline the analysis required. The kernel of the Hermitian operator is

p (\tilde{r} - \tilde{r}', \tilde{z}, \tilde{z}') = {(\frac{s - \tilde{z}}{s - \tilde{z}'})}^{2} \sum_{m - 1}^{M} δ (\tilde{r} - \tilde{r}' + r_{m} \frac{\tilde{z} - \tilde{z}'}{s}),

leading to the object-space eigenequation

\int_{{\tilde{z}}_{1}}^{{\tilde{z}}_{2}} d \tilde{z}' {(\frac{s - \tilde{z}}{s - \tilde{z}'})}^{2} U_{j} (ρ, \tilde{z}') \sum_{m = 1}^{M} exp (2 π i ρ \cdot r_{m} (\tilde{z} - \tilde{z}') \frac{1}{s}) = μ_{ρ, j} U_{j} (ρ, \tilde{z})

where it has been assumed the singular functions can be written as

u_{j, ρ} (\tilde{r}, \tilde{z}) = U_{j} (ρ, \tilde{z}) exp (2 π i ρ \cdot \tilde{r}) .

Introducing Û(ρ,z̃)(s − z̃)² = U(ρ,z̃) we find the integral equation

{\hat{U}}_{j} (ρ, \tilde{z}) = \int_{{\tilde{z}}_{1}}^{{\tilde{z}}_{2}} d \tilde{z}' \hat{U} (ρ, \tilde{z}') \sum_{m = 1}^{M} exp (2 π i \frac{\tilde{z}}{s} ρ \cdot r_{m}) exp (- 2 π i \frac{\tilde{z}'}{s} ρ \cdot r_{m}) .

This integral equation has a degenerate kernel, and such equations can be solved using standard methods [19].

Equation (19) is the kernel of the object-space Hermitian operator, but can also be seen the impulse-response function of the system. A delta function at (r̃′,z̃′) in the object will cause a distribution p(r̃ − r̃′,z̃,z̃′) in the reconstruction. The impulse response function is a sum of M lines whose intensity varies with z̃, and that all meet at r̃ = r̃′, z̃ = z̃′.

3. Numerical results

In section 2, the singular functions have been found for the simplified case of infinite object and image. These results can be illustrated numerically. Under certain circumstances, the SVD technique is not very time-consuming and can be run on a normal desk-top computer. The calculation times are on the order of minutes or tens of minutes. The most time-consuming part is actually the simulation of the object (see section 3.2). The technique for numerical SVD calculations is very similar to Burvall et al. [14].

3.1. Numerical SVD

The geometry is the same as in Fig. 1. We have assumed a distance of 30 cm from source to image plane, and that the object is located between 10 and 1 cm from the image plane, i.e. s = 0.3 m, z ₁ = −0.1 m, and z ₂ = −0.01 m. Furthermore, it is assumed that all the sources are placed along a line or arc rather than in a 2D array. While it would be interesting to investigate other arrangements, this option simplifies the analysis significantly and is often employed in tomosynthesis. Assuming the source distribution is along the x axis, we have r _m = (ξ_m, 0). The outermost source positions are at χ_m = ±0.1 m, and the source positions evenly distributed in between. The examples below are for 11 source positions, i.e., M = 11.

The assumption that r _m = (χ_m, 0) is used to simplify the analysis by reducing large parts of it from 3D to 2D. A closer look at Eq. (18) reveals that it depends on ρ only as ρ · (r _m − r _m′) = ρ_x(χ_m − χ _m′), so the elements of $K$ depend only on ρ_x and not on ρ_y. Comparison to Eq. (16) shows that V _j(ρ) = V _j(ρ_x) or μ_ρ,j = μ_{ρ_x,j} will not depend on ρ_y either.

Once the image-space singular functions are known, the object-space singular functions U_ρ,j(r̃,z̃) can be found from the projection

[H^{†} v_{ρ, j}] = \sqrt{μ_{ρ, j}} u_{ρ, j} (\tilde{r}, \tilde{z}) .

This equation is also known as part of the shifted-eigenvalue equations [20]. Insertion of Eqs. (9), (8) and (15) into Eq. (20) yields

u_{ρ, j} (\tilde{r}, \tilde{z}) = exp (2 π i ρ \cdot \tilde{r}) U_{j} (ρ, \tilde{z})

where

U_{j} (ρ, \tilde{z}) = \frac{1}{\sqrt{μ_{ρ, j}}} {(\frac{s - \tilde{z}}{s})}^{2} \sum_{m = 1}^{M} {[V_{j} (ρ)]}_{m} exp (2 π i \frac{\tilde{z}}{s} ρ \cdot r_{m}) .

So U_j(ρ,z̃) = U_j(ρ_x,z̃) only needs to be calculated for ρ = (ρ_x, 0) and then the same 2D matrix can be used for all ρ_y. Figure 3 shows the eigenvalues μ_{ρ_x,j} while Fig. 4 contains the absolute values of some of the singular functions U_j(ρ_x,z̃). We note that the eigenvalues are approximately zero for low frequencies, and that this band of zero eigenvalues grows broader with increasing j. The eigenvalues are a combination of discrete and continuous spectra: discrete in j which handles the axial dimension, but continuous in ρ_x which handles the transverse dimensions.

Fig. 3 Eigenvalues μ_{ρ_x,j} for j ranging from 1 to 11.

Download Full Size | PDF

Fig. 4 Absolute values of the singular functions (a) U ₁(ρ_x,z̃), (b) U ₂(ρ_x,z̃), (c) U ₆(ρ_x,z̃), and (d) U ₁₁(ρ_x,z̃).

Download Full Size | PDF

Unlike the 3D imaging case [14], where the singular functions are real, these functions are complex. In Fig. 4 their absolute value is shown. Interpreting the results intuitively is fairly difficult, as there are no obvious image planes and the coordinates are a combination of spatial (z̃) and frequency (ρ). It is easier to interpret the expansion of an example object into these eigenfunctions, i.e., its measurement component. This requires a relevant sample object.

3.2. Sample object

We illustrate the method using an example of at least some relevance to the future task, namely tomosynthesis of breast tissue. The object consists of two parts: a designer nodule and a clustered lumpy background (CLB). This approach is chosen since the expected advantage of tomosynthesis is the ability to detect low-contrast masses, such as the nodule, that are masked by overlapping tissue, such as the clustered lumpy background. The designer nodule has a projection onto the image plane given by given by [21]

g_{m} (r_{d}) = A rect (\frac{2 | r_{d} |}{R}) {(1 - \frac{{| r_{d} |}^{2}}{R^{2}})}^{v}

where A is the amplitude of the nodule, R its radius, and v a real and non-negative shape factor that determines the sharpness of the edges. We have chosen R = 3mm, A = 0.03, and v = 1.5. This value of v is reported to give the best fit to average tumor profiles [21]. The 3D object that gives this projection is [22]

f (r) = A \frac{3}{4 R} rect (\frac{2 | r |}{R}) (1 - \frac{{| r |}^{2}}{R^{2}}),

as can be shown by finding its projection for a far-away source, i.e., by integrating Eq. (22) with respect to z.

The CLB simulates the 2D projection of the background by distributing lumps according to a statistical model [23]. The projection is given by

g (r_{d}) = \sum_{k = 1}^{K} \sum_{n = 1}^{N_{k}} b (\frac{1}{a_{k n}} (r_{d} - r_{k} - r_{k n}), R_{θ_{k} n})

where K is the numbers of clusters, each centered on r _k, and N_k is the number of blobs or lumps in each cluster, their positions within the cluster given by r _kn. The function b describes the shape of each blob, a_kn its strength, and R _θ the rotation matrix of the random angle θ_kn. The theory of this 2D background is thoroughly described by Bochud et al. [23]. However, no proper 3D model exists and we need to model the 3D object rather than its projection. Since the intention in the current paper is only to provide an illustration, we have extended the model in the simplest possible way. The lumps are made three-dimensional, and placed in a volume by the 3D vectors r _k and r _kn. They were also rotated in 3D rather than in 2D. Apart from that, we followed the classic CLB as introduced by Bochud [23]. The number K was a Poisson random variable with expectation value 150 and N_k a Poisson random variable with expectation value 20. The scale a_kn was set to 1 for all blobs, and the position of the cluster r _k a uniform random variable within the object area which is 40 × 40 × 90mm. The position within the cluster r _kn is a Gaussian with standard deviation 3.6mm. The function b is an asymmetrical exponential blob b(r,R _{θ_kⁿ}) = exp(−α|R _{θ_kⁿ} r|^β/L(R _{θ_kⁿ} r)) as given by Bochud [23], of length L_x = 15mm and width L_y = L_z = 6mm. Parameters of the blob function are α = 2.1 and β = 0.5. The blob is rotated an angle θ_kn = θ_k (uniform random variable on [0,2π]) in the xy-plane, and an angle ϕ_kn = ϕ_k (uniform random variable on [0, π]) in the xz-plane.

An example object f ̃(r̃,z̃) is shown in Fig. 5, in the distorted coordinates r̃ and z̃. It is shown in two slices: one through the middle of the designer nodule in the x̃ỹ plane, and one through the middle of the designer nodule in the x̃z̃ plane. The object is not a proper 3D breast model, since it is not connected enough for a breast CT slice [24], but it will give reasonable results once projected and will serve as an illustration.

Fig. 5 Example object f ̃(r̃, z̃), including a designer nodule and a 3D clustered lumpy background. (a) Section through the object at z̃ = −60mm, or z = −50mm. (b) Section through the object at ỹ = ỹ = 0. Both sections go through the center of the designer nodule.

Download Full Size | PDF

3.3. Numerical results

Once an object has been generated, and the singular functions of the system found, they can be combined to generate the measurement component. It is an expansion of the object into the object-space singular functions according to

{\tilde{f}}_{meas} (\tilde{r}, \tilde{z}) = \iint_{- \infty}^{\infty} d^{2} ρ \sum_{j = 1}^{M} A_{j} (ρ) U_{j} (ρ, \tilde{z}) exp 2 π i ρ \cdot \tilde{r},

where the integral is numerically found as a Fast Fourier Transform. The expansion coefficients A_j(ρ) are found as

A_{j} (ρ) = \int_{{\tilde{z}}_{1}}^{{\tilde{z}}_{2}} d \tilde{z} \tilde{F} (ρ, \tilde{z}) U_{j}^{*} (ρ, \tilde{z}) \frac{s^{4}}{{(s - \tilde{z})}^{4}}

where F̃(ρ,z̃) is the 2D Fourier transform of the object f̃(r̃,z̃) with respect to r̃. Two things should be noted, as they differ from the expression used in Ref. [14]. First, the complex conjugate in Eq. 23 was omitted in [14] since the singular functions were real, but must be included here. Second, the weighting function w(x̃,z̃) = s ⁴/(s − z̃)⁴ follows from the earlier change of variables. Slices through the measurement component are shown in Fig. 6. The slices are taken at the same positions as those in Fig. 5. In part (b), the low axial resolution is showing up as elongation of both nodule and background.

Fig. 6 Measurement component of the object in Fig. 5. (a) Section through the measurement component at z̃ = −60mm, or z = −50mm. (b) Section through the measurement component at ỹ = y = 0. Both sections go through the center of the designer nodule.

Download Full Size | PDF

The zero-angle projection of the object is shown in Fig. 7(a), along with the same projection of its measurement component in part (b). The images are generated by direct numerical propagation of the quantities involved, using a sampling of 400 points in the transverse directions and 300 points in the longitudinal direction. This serves as a test of both the analysis and the numerical accuracy: the image of an object should be identical to the image of its measurement component. As can be seen in Fig. 8, the difference between them is below 1% except at the edges were edge effects could cause it to rise to around 10%.

Fig. 7 (a) Image or projection of the measurement component in Fig. 6. (b) Image or projection of the object in Fig. 5. Both are calculated through direct numerical propagation of the object and measurement component respectively, for the source placed on the z-axis.

Download Full Size | PDF

Fig. 8 Difference between the projection of the measurement component (see Fig. 7(a)) and the projection of the object (see Fig. 7(b)). The maximum absolute value is around 3 × 10⁻⁶ compared to the maximum of the projections, which is around 2 × 10⁻⁴. This gives an error of roughly 1%.

Download Full Size | PDF

An interesting point is how the slice through the measurement component is a hybrid between the projection and the similar object slice. This illustrates the very principle of tomosynthesis. Unlike tomography, which provides a full reconstruction of every plane of the object, tomosynthesis provides only partial 3D information. The longitudinal resolution will always be low, and contributions from close-by transverse planes will be present in every slice of the measurement component.

4. Discussion and future work

The measurement component of an object provides information on how well this object is transmitted through the imaging system. Comparing the measurement component in Fig. 6 to its object in Fig. 5, we can see that the information received in the image plane is not complete. This illustrates the fundamental limit of the acquisition process, and serves as an estimate of the best reconstruction that could be achieved. We can also see that the measurement component looks a lot like the object, but that it also borrows traits from the projection in Fig. 7. This reflects the small angle used in tomosynthesis and its effect on the results. The reconstruction, while given in 3D, is in many ways an improved projection or a projection with some separation of information into different layers. It is not the full 3D reconstruction that could be achieved from e.g. tomography.

It is also becoming evident that there is a lack of three-dimensional models of the breast tissue. Statistical models for the normal breast, such as the clustered lumpy background [23], are made for 2D and not for 3D. Since mammography has been the dominating investigation method, it has simply been sufficient to work only on the 2D projection through the breast. Extensions to 3D are either intuitive and lack the statistical back-up required for a systematic study, like the example presented in this paper, or advanced models designed to produce one really good phantom and not the large number of samples required for a statistical study. There is need for a model of the breast tissue that has the same statistical properties as experimental images in 3D but also in its 2D projections, that looks similar to experimental images in 3D but also in its 2D projections, and that allows for generation of numerous realizations. Strong candidates at the moment are a more sophisticated extension of the clustered lumpy background [23], the threshold model by Abbey and Boone [25], or the breast phantoms by Zhang et al. [26]. The example object in Fig. 5 is sufficient as an illustration of the SVD method, but is not interconnected enough for a real slice through breast tissue.

The next step would be to use the SVD of a more realistic breast tissue model in an image-quality study, including the effect of noise. The SVD model is particularly suitable for this since the statistics of the object background are relatively easily transferred to the image through the use of characteristic functionals [11, Sec. 8.5.3]. The SVD singular functions can also be used as channels for a channelized observer [27]. This allows for many comparisons, for example how or if the detection seems to improve with tomosynthesis compared to mammography, or how the SVD method compares to other reconstruction algorithms.

5. Conclusions

The mixed continuous and discrete operators of a tomosynthesis system have been derived from the Boltzmann equation, and used to develop an analytical singular-value decomposition of a tomosynthesis system. In combination with numerical calculations, it provides the measurement component of any object, where object refers to the breast tissue. The measurement component of an example object has been presented. Although the example object is not truly realistic, interesting information on the tomosynthesis process and its effects on reconstruction have been found. Future prospects involve more realistic 3D breast tissue models, and an image-quality study involving a channelized Hotelling observer.

Acknowledgments

C. Dainty acknowledges support from Science Foundation Ireland Grant No 07/IN.1/1906, H. Barrett from the Center for Gamma-Ray Imaging Grant 5 P41EB002035-11 and the SPECT Imaging and Parallel Computing Grant 5 R37 EB000803-19, and A. Burvall from the Swedish Research Council.

References and links

1. J. T. Dobbins III and D. J. Godfrey, “Digital x-ray tomosynthesis: current state of the art and clinical potential,” Phys. Med. Biol. 48, R65–R106 (2003). [PubMed]

2. J. T. Dobbins III, “Tomosynthesis: at translational crossroads,” Med. Phys. 36, 1956–1967 (2009). [PubMed]

3. L. T. Niklason, “Digital tomosynthesis in breast imaging,” Radiology 1997, 399–406 (1997).

4. G. Gennaro, A. Toledano, C. di Maggio, E. Baldan, E. Bezzon, M. La Grassa, L. Pescarini, I. Polico, A. Proietti, A. Toffoli, and P. C. Muzzio, “Digital breast tomosynthesis versus digital mammography: a clinical performance study,” Eur. Radiol. (2009). [PubMed]

5. I. Andersson, D. M. Ikeda, S. Zackrisson, M. Ruschin, T. Svahn, P. Timberg, and A. Tingberg, “Breast tomosynthesis and digital mammography: a comparison of breast cancer visibility and BIRADS classification in a population of cancers with subtle mammographic findings,” Eur. Radiol. 18, 2817–2825 (2008). [PubMed]

6. W. F. Good, G. S. Abrams, V. J. Catullo, D. M. Chough, M. A. Ganott, C. M. Hakim, and D. Gur, “Digital breast tomosynthesis: a pilot observer study,” Am. J. Radiology 190, 865–869 (2008).

7. S. P. Poplack, T. D. Tosteson, C. A. Kogel, and H. M. Nagy, “Digital breast tomosyntheis: Initial experiance in 98 women with abnormal digital screening mammography,” Am. J. Radiology 189, 616–623 (2007).

8. A. S. Chawla, J. Y. Lo, J. A. Baker, and E. Samei, “Optimized image acquisition for breast tomosynthesis in projection and reconstruction space,” Med. Phys. 36, 4859–4869 (2009). [PubMed]

9. T. Wu, R. H. Moore, E. A. Rafferty, and D. B. Kopans, “A comparison of reconstruction algorithms for breast tomosynthesis,” Med. Phys. 31, 2636–2647 (2004). [PubMed]

10. Y. Zhang, H.-P. Chan, B. Sahiner, J. Wei, M. M. Goodsitt, L. M. Hadjiiiski, J. Ge, and C. Zhou, “A comparative study of limited-angle cone-beam reconstruction methods for breast tomosynthesis,” Med. Phys. 33, 3781–3795 (2006). [PubMed]

11. H. H. Barrett and K. J. Myers, Foundations of Image Science (John Wiley, Hoboken, New Jersey, 2004).

12. M. Bertero and P. Boccacci, Inverse Problems in Imaging (Institute of Physics Publishing, Bristol, UK, 1998)

13. H. H. Barrett, J. N. Aarsvold, and T. J. Roney, “Null functions and eigenfunctions: tools for the analysis of imaging systems,” Lect. Notes. Comput. Sci. 11, 211–226 (1991).

14. A. Burvall, H. H. Barrett, C. Dainty, and K. J. Myers, “Singular-value decomposition for through-focus imaging systems,” J. Opt. Soc. Am. A 23, 2440–2448 (2006).

15. J. Yao and H. H. Barrett, “Predicting human performance by a channelized Hotelling observer,” Proc. SPIE 1768, 161–168 (1992).

16. H. H. Barrett, J. Yao, J. P. Rolland, and K. J. Myers, “Model observers for assessment of image quality,” Proc. Natl. Acad. Sci. USA 90, 9758–9765 (1993). [PubMed]

17. M. Y. Chiu, H. H. Barrett, R. G. Simpson, C. Chou, J. W. Arendt, and G. R. Gindi, “Three-dimensional radio-graphic imaging with a restricted view angle,” J. Opt. Soc. Am. 69, 1323–1333 (1979).

18. R. Pierri, A. Liseno, F. Soldovieri, and R. Solimene, “In-depth resolution for a strip source in the Fresnel zone,” J. Opt. Soc. Am. A 18, 352–359 (2001).

19. A. D. Polianin and A. V. Manzhirov, Handbook of integral equations (CRC Press, Florida, 1998) chapter 11.2.

20. C. Lanczos, Linear differential operators (Van Nostrand, London, 1961).

21. A. E. Burgess, “Visual signal detection with two-component noise: low-pass spectrum effects,” J. Opt. Soc. Am. A 16, 694–704 (1999).

22. I. Reiser and R. M. Nishikawa, “Task-based assessment of breast tomosynthesis: Effects of acquisition parameters and quantum noise,” Med. Phys. 37, 1591–1600 (2010). [PubMed]

23. F. O. Bochud, C. K. Abbey, and M. P. Eckstein, “Statistical texture synthesis of mammographic images with clustered lumpy backgrounds,” Opt. Express 4, 33–43 (1998).

24. K. G. Metheany, C. K. Abbey, N. Packard, and J. M. Boone, “Characterizing anatomical variability in breast CT images,” Med. Phys. 35, 4685–4694 (2008). [PubMed]

25. C. K. Abbey and J. M. Boone, “An ideal observer for a model of x-ray imaging in breast parenchymal tissue,” (E.A. Krupinski, Ed.): IWDM 2008, LNCS5116, 393–400 (2008).

26. C. Zhang, P. R. Bakic, and A. D. A. Maidment, “Development of an anthropomorphic breast software phantom based on region growing algorithm,” Proc. SPIE 6918, 69180V (2008).

27. S. Park, J. M. Witten, and K. J. Myers, “Singular vectors of a linear imaging system as efficient channels for the bayesian ideal observer,” IEEE Trans. Med. Imaging 28, 657–668 (2009). [PubMed]

Singular-value decomposition of a tomosynthesis system

Abstract

1. Introduction

2. SVD analysis

2.1. Infinite object and image

3. Numerical results

3.1. Numerical SVD

3.2. Sample object

3.3. Numerical results

4. Discussion and future work

5. Conclusions

Acknowledgments

References and links

Cited By

Figures (8)

Equations (39)

Optics Express