Imaging of Fiber-Like Structures in Digital Breast Tomosynthesis

Data for both physical phantom and patient studies was acquired with a Hologic Selenia Dimensions DBT scanner, in which 15 projections are obtained over a scanning arc which is approximately 15 degrees. With the actual scanner, both X-ray source and flat-panel detector move with each projection. The geometric configuration data are recorded with each scan and are stored in the form of projection matrices. These matrices map 3D spatial coordinates of the volume to 2D detector coordinates for each projection view Li, Da, and Liu (2010).

To explain the DBT data model and image reconstruction algorithms, we employ an approximate simplified geometric model for the DBT scan shown in Fig. 1. In this model, the detector is fixed; the X-ray source moves on a circular arc of 15 degrees, acquiring 15 projections in 1 degree increments; the center of the scanning arc lies on the detector; and the detector pixel size is 0.14 mm $\times$ 0.14 mm. The various geometric parameters of the scan are also shown in the schematic in Fig. 1.

While the data model and image reconstruction theory is explained with the approximate stationary-detector model, the image reconstruction results are obtained with the precise scan geometry provided by the projection matrices.

The X-ray transmission data are modeled with the Beer-Lambert law

where $I_{0}$ is the incident X-ray fluence; $I(\xi,u,v)$ is the transmitted fluence for X-ray source angle $\xi$ and detector pixel location $(u,v)$. The function $g(\xi,u,v)$ represents the line integrals over the object function $f(\vec{r})$

where $f(\vec{r})$ represents an energy-averaged map of the linear X-ray attenuation coefficient; the X-ray source position is $\vec{r}_{0}(\xi)$; and the unit vector $\hat{\theta}(\xi,u,v)$ indicates the direction of a ray originating from $\vec{r}_{0}(\xi)$ and incident on the detector at $(u,v)$. Using the coordinate system shown in Fig. 1, the source position can be parameterized by

and the detector pixel locations by

where $C=70$ cm, the radius of the source trajectory, and $L=28.7$ cm, the length of the detector in the $x$ direction. The unit vector $\hat{\theta}(\xi,u,v)$ can then be written

The $u$ coordinates are shifted by $L/2$ with respect to the volume coordinates because the origin of the detector is taken to be at a corner, while the origin of the volume Cartesian coordinates coincides with the middle of the detector edge adjacent to the chest wall.

The noisy simulation data, which is generated for the results in Sec. III, use the Beer-Lambert model in Eq. (1) as a mean value for a Poisson distribution, where the noise level is specified by selecting the total number of photons, $I_{0}$, incident at each detector bin. No correlation between detector bins is assumed.

For image reconstruction, the DBT transmission data are processed by inverting Eq. (1)

to obtain the projection data. Image reconstruction algorithms are based on the mathematical inversion of the X-ray projection model in Eq. (2).

For the results presented in Sec. III, four image reconstruction algorithms are compared. The first two are of the form of filtered back-projection (FBP). For background on application of FBP in CT, the reader is referred to Ref. Kak and Slaney, 1988, and for implementation of FBP specific to DBT, Ref. Mertelmeier, 2017. We also consider two common forms of penalized, least-squares optimization that we have already investigated in the context of signal detection in DBT Sanchez, Sidky, and Pan (2015); Rose et al. (2017).

For the FBP implementations, the general form of filtering considered involves only filtering along the detector $u$-direction with the one-dimensional Fourier transform (FT) and inverse Fourier transform (IFT)

where integration over $u^{\prime}$ performs the FT; $H(\nu)$ is the filter function; integration over $\nu$ performs the IFT; and $h(\xi,u,v)$ is the filtered projection data.

The reconstructed volume is obtained from $h(\xi,u,v)$ by back-projection,

where $u(\xi,\vec{r})$ and $v(\xi,\vec{r})$ are the coordinates on the detector pointed to by a ray originating at source position $\vec{r}_{0}(\xi)$ and passing through $\vec{r}=(x,y,z)$, i.e.

FBP - The filter for the first FBP algorithm is a ramp filter with a Hanning apodizing window

where $\nu_{H}$ is the cutoff frequency of the filter. The cutoff frequency $\nu_{H}$ is specified as a fraction of the detector’s Nyquist frequency $\nu_{max}$

where $\Delta u$ is the detector bin size and $\nu_{c}$ is the fractional cutoff frequency. For the studies in Sec. III, we investigate the impact of $\nu_{c}$ on conspicuity of fiber-like signals.

mFBP - One strategy for designing analytic image reconstruction for tomosynthesis involves combining back-projection BP with ramp-filtered FBP. Such an admixture has been proposed for reduction of metal artifacts Gomi, Hirano, and Umeda (2009). A similar approach is to modify the FBP filter in such a way that its value at low spatial frequency is boosted with respect to the ramp Orman, Mertelmeier, and Haere (2006). For mFBP, short for modified FBP, we design a filter based on the latter idea. Our particular choice of filter function employs a shifted quadratic at low spatial frequency

where $\delta$ is specified in terms of a fraction of the detector’s Nyquist frequency. The filter $H_{\text{mFBP}}$ used for mFBP is the product of the boosted filter $H_{\delta}(\nu)$ and the Hanning apodization window $H_{\text{Hanning}}(\nu)$ for $\nu_{c}=1.0$. For mFBP, $\delta$ is the only regularization parameter and is taken to be in the range $[0,\pi^{-1}]$.

At $\delta=0$, $H_{\text{mFBP}}(v)$ is identical to $H_{\text{FBP}}(v)$ for $\nu_{c}=1$. The upper bound on the range $\delta=\pi^{-1}$ is determined by the condition $\frac{\mathrm{d}^{2}H_{\text{mFBP}}}{\mathrm{d}\nu^{2}}(0)=0$. This condition yields a flat response at low frequencies (the first derivative at $\nu=0$ is zero because the filter function is symmetric). At $\delta=\pi^{-1}$, $H_{\text{mFBP}}(v)$ becomes a low-pass filter and the corresponding mFBP reconstruction is approximately back-projection.

We point out that there is other work on deriving filters for FBP applied to tomosynthesis Stevens, Fahrig, and Pelc (2000); Kunze et al. (2007); Nielsen et al. (2012). Also, additional parameters such as a slice thickness filter or Hanning filter width can be beneficial for controlling the DBT image quality. Orman, Mertelmeier, and Haere (2006); Mertelmeier (2017) Our particular design results from empirical subjective assessment on volumes reconstructed with the Hologic geometry and the aim of having a single-parameter image reconstruction algorithm.

LSQI - We also investigate whether orientation-dependence of fiber-like signal conspicuity can be mitigated using iterative reconstruction algorithms. Such algorithms employ a discrete model, and for the present work we employ a voxel expansion to arrive at a discretized form of Eq. (2)

where $f$ is a vector of $m$ voxel coefficients; $g$ is a vector of the $n$ line-integral measurements; and $X$ is the system matrix generated by computing the projection of each of the voxel expansion elements. The dimensions of the voxels used here measure 0.14$\times$0.14 mm² in-plane and 1.3 mm in depth. In using $g$ and $f$ without parenthetical arguments, we are referring to the discrete digital data and voxel coefficients, respectively. With parenthetical arguments, $g(\xi,u,v)$ and $f(\vec{r})$, we are referring to the continuous data function and volume, respectively.

We consider two forms of least-squares optimization. The first form of least-squares employs identity Tikhonov regularization (LSQI)

where $\lambda$ is a normalized regularization parameter and $\|X\|$ is the maximum singular value of the system matrix $X$. For more implementation details on LSQI for DBT, consult Ref. Rose et al., 2017.

In this particular formulation of LSQI, the parameter $\lambda$ is dimensionless and the same value of $\lambda$ will incur the same level of regularization no matter what length units are used in specifying $X$. Some sense of the impact of $\lambda$ can be deduced by deriving the small and large $\lambda$-value limits of the LSQI analytic solution.

Setting the gradient of the objective function in Eq. (5) to zero, and solving for $f$ yields a direct linear expression relating $f_{\text{recon}}$ to the data $g$

This expression is useful for interpreting the effect of the regularization parameter $\lambda$. The transpose of the projection matrix $X^{T}$ represents a form of discrete back-projection. The small $\lambda$ limit tends to back-projection filtration (BPF)

where the discrete back-projection, $X^{T}$, is followed by filtration in the form $(X^{T}X)^{-1}$. Scaled appropriately, the large $\lambda$ limit tends to back-projection Rose et al. (2017)

In practice, $f_{\text{recon}}$ is obtained by iterative solution of Eq. (5) using the conjugate gradients algorithm.

LSQD - The second form of least-squares presented is more commonly used in tomographic applications than LSQI; namely, we consider least-squares with a quadratic roughness penalty Fessler (1994) (LSQD)

where $\nabla$ is a finite differencing matrix approximating the gradient operator. The scale factor $S=\|X\|/\|\nabla\|$ serves to make the regularization parameter $\lambda$ dimensionless.

The direct expression for the LSQD volume is

The small $\lambda$ limit of this expression is the same as that of LSQI, a form of back-projection filtration. The scaled large $\lambda$ limit, however, differs

The expression $(\nabla^{T}\nabla)^{-1}$ is a form of blurring, and as a result the large $\lambda$ limit of LSQD is blurred back-projection.

Pyramidal image volume - The Hologic system reconstructs images onto a pyramidal volume to correct for the effects of magnification due to the extremely limited angular range. In terms of the continuous reconstructed volume, slices are displayed in the transformed coordinate system $(x^{\prime},y^{\prime},z^{\prime})$, where the transformed coordinates are related to the original Cartesian coordinates by

The reconstructed volume is represented in this new coordinate system by the function

Note that the transformed coordinates match the untransformed coordinates at the plane of the detector, i.e. $z=0$. As $z$ increases, $x^{\prime}$ and $y^{\prime}$ adjust to counteract the magnification of the corresponding $z$-plane. Results for all four investigated algorithms are presented in the transformed coordinates.

Background flattening - Raw reconstructions of the breast exhibit a low-frequency variation when reconstructing at larger regularization strengths. The effect is particularly strong near the skin line. This gray-level variation needs to be removed in order to display low-contrast tissue structures over extended areas of the breast. In the present work, we fit a low-order polynomial to the 2D slice images then divide the raw 2D slice image by this fitted polynomial.

In order to provide background on fiber-like signal imaging in DBT and to motivate the relative intensity metric, we show tomosynthesis images of a rod aligned parallel and perpendicular to the direction of X-ray source travel (DXST) in Fig. 2. Shown are images for FBP using the ramp filtering

and no filtering

which corresponds to back-projection (BP). Because of the extremely limited angular range of the scan trajectory, the reconstructed slice images resemble closely the central filtered projection prior to back-projection.

In the FBP image of Fig. 2, the ramp filter causes two classic features of objects reconstructed in a tomosynthesis settingReiser and Glick (2014); Sechopoulos (2013): the ramp filter amplifies high frequency content in the $x$-direction causing the edges perpendicular to the DXST to enhance, and it attenuates low frequency content in the $x$-direction causing gray-level contrast to drop toward the middle of the signal. The impact of these FBP features differ greatly between the images of the two rods. For the rod parallel to the DXST, the edge-enhancement is minimal because the short edges are perpendicular to the DXST; furthermore, the gray-level contrast is low relative to the actual rod because the long-axis is parallel to the DXST. For the rod perpendicular to the DXST, the situation is reversed and the rod enhances well. These basic imaging properties are a direct result of the ramp filter, and this is seen in comparing the DBT slice images with the ramp filtered central projection in Fig. 2.

The BP image shown in Fig. 2 does not have either the edge-enhancement nor the drop in contrast in the signal interior. The central slice of the BP image appears to depict both rods faithfully. In comparing the FBP to BP images, the former clearly have a greater variation with in-plane orientation with respect to the DXST. In terms of conspicuity, BP will likely yield uniform results with respect to orientation. For FBP, conspicuity may vary greatly with orientation.

While the in-focus-plane images of Fig. 2 would seem to favor BP, there is, however, a penalty for the absence of filtering; namely the depth-blur of BP image reconstruction is significantly greater than that of FBP. The depth blur of reconstructed rod signals of different lengths is shown in Fig. 3 in a non-standard viewing plane. The images are shown for planes aligned along the DXST and perpendicular to the detector. The BP image strongly displays the classic size-dependent depth blur of DBT; the BP depth blur has a central diamond shaped region of nearly uniform gray values, and outside of this region the intensity slowly dissipates. The geometry of the central diamond region is determined by the $x$-extent of the rod and the scanning angular range, which for the present configuration is 15^∘. For FBP, the reconstructed rod images in the depth plane have quite different appearances from those of BP. The ends of the rod are highlighted, but for the longer rods the “BP diamond” is suppressed including the mid-line where the actual rod is located.

Results from the investigated image reconstruction algorithms can be interpreted as some combination of the two extreme cases of FBP and BP. There is a trade-off for fiber-like signals. For FBP-like algorithms, there is strong orientation-dependence of the conspicuity in the focus plane. For BP-like algorithms, there is pronounced depth blur for fibers parallel to the DXST. To quantify this trade-off we employ a relative intensity metric defined through noiseless simulation of rod-signals.

Zoomed in reconstructions of a subset of the rods are shown in Fig. 5 along with the orientation-dependent relative intensity metric. The shown images are in-plane and the relative intensity metric is abbreviated as $i(\phi)$ with $z=0$ and $L=1.9$ mm. $i(\phi)$ captures the trends in visual rod intensity relative to that of $\phi=90^{\circ}$ for the various algorithms and regularization strengths. Overall for each of the image reconstruction algorithms the deviation of $i(\phi)$ from 1.0 is greatest for low regularization strength and small $\phi$. Increasing regularization reduces the orientation-dependent rod intensity as the low $\phi$ rods show values of $i(\phi)$ approaching 1.0.

The detailed behavior of each algorithm differs as each uses a different form of regularization. The FBP Hanning filter regularizes by smoothing the image in the $x$-direction, and orientation-dependent intensity is mitigated with the smoothing filter width, defined by $\nu_{c}$, comparable to the length of the rod signal. While the intensity of $0^{\circ}$ and $90^{\circ}$ rods appears similar at large regularization (low $\nu_{c}$), the distortion of the rod has strong $\phi$ dependence because of the directionality of the filter. The mFBP filter parameter $\delta$ interpolates between ramp-FBP and BP, and the rod signal intensity and shape become more uniform with increasing $\delta$. The LSQI algorithm parameter $\lambda$ ranges between similar extremes as mFBP and accordingly the resulting rod images are similar to mFBP. For LSQD, the low $\lambda$ limit is similar to that of LSQI, but the large $\lambda$ limit involves smoothing of a different form than the other three algorithms, which is reflected in the different appearance of the rod images.

The next set of simulation results illustrate how the intensity variation of fiber-like signals reconstructed from noiseless data translate to variation in signal conspicuity in the presence of noise. The simulated data are generated from projections with the complete test phantom illustrated in Fig. 4 that contains a uniform background block and 10 rod signals oriented in angular increments of $10^{\circ}$ over a range of $0^{\circ}$ to $90^{\circ}$ with respect to the $x$-axis. Noise is modeled as an independent Poisson process in the transmission data with a mean incident intensity of $I_{0}=1.5\times 10^{4}$ photons per ray.

The central slice from each reconstructed volume is shown in Fig. 6. Assessing the shown images for rod conspicuity yields an overview of in-plane orientation-dependence and the impact of the various algorithms and forms of regularization. Starting with FBP in the upper left set of images, there is a clear orientation-dependence of the rod conspicuity for $\nu_{c}=1.0$, the least regularized case. The rod conspicuity increases with increasing angle to the $x$-axis, the DXST. Increasing regularization, by decreasing $\nu_{c}$, visually seems to improve the conspicuity of the low-angle rods.

For each of the other three algorithms, mFBP, LSQI, and LSQD, at an overview level, we observe similar trends to that of FBP for rod conspicuity. For the least regularized phantom image, there is a strong orientation dependence with rod angle. Increasing regularization strength improves the conspicuity of the rod-signals nearly parallel to the DXST. These rod-conspicuity results are fairly consistent despite the obvious differences in other image quality characteristics, such as noise texture or blur.

To demonstrate the in-plane anisotropy of fiber-like signals, we show FBP reconstructions ($\nu_{c}=1.0$) of the ACR-DM phantom at two orientations in Fig. 7. In Fig. 7 the DXST is from left to right, and in the top row the ACR-DM phantom is oriented at approximately a 40^∘ angle from the DXST. Orienting the phantom in this way causes the odd-numbered fibers to line up nearly perpendicular to the DXST. The top-left panel image shows the first three fibers, and the striking feature is that the second fiber is invisible while fibers one and three are clearly seen. Even more striking is that the fibers are numbered in order of decreasing thickness, yet, fiber number three is seen and thicker fiber number two is not. The bottom-row illustrates the imaging properties of the same phantom when tilted 10^∘ from the top-row orientation. This change in ACR-DM phantom orientation causes the second fiber to appear in the bottom-left panel image. The top-row image clearly displays how fiber-like signals enhance when oriented perpendicular to the DXST while those parallel do not; the bottom-row image hints at the strong dependence of fiber-like signal conspicuity for orientations close to the DXST.

We expand upon the results shown in Sec. III.3 for the ACR-DM phantom tilted at $40^{\circ}$ so that the fiber objects line up at $5^{\circ}$ and $95^{\circ}$ with respect to the DXST. The ACR-DM phantom data is reconstructed by the four considered algorithms at matched resolution with FBP. The regularization equivalence among the four reconstruction algorithms is established by matching the artifact spread function (ASF) Wu et al. (2004); Hu, Zhao, and Zhao (2008) of a simulated small point-like object: a sphere of 28.0 $\mu$m diameter and attenuation coefficient of 9.29 cm^-1, modeling a microcalcification.

In Fig. 8, we show ACR-DM images for the 40^∘ orientation so that fiber signals at 5^∘ and 95^∘ are presented. The 5^∘ fiber is difficult to see; it is visible in three of the mFBP images, two of the LSQI images, and none of the LSQD and FBP images . We note that there are differences in the noise texture among the images from the various algorithms. In particular, that of LSQD is markedly different than that of the other three, and the differences in noise texture can impact the fiber conspicuity.

Fig. 9 is identical with Fig. 8 except that the shown images are from a slice 7.7 mm below the focus plane for the ACR-DM fibers. The shown plane is well below the fibers. The main point of viewing this out-of-focus plane is to observe the extent to which fiber signals bleed through, yielding a sense of the object-dependent blur. Most notable is the fact that the 5^∘ fiber is still visible for the more regularized images of mFBP and LSQI. Thus, it appears that the improved conspicuity of the $5^{\circ}$ fiber for the regularized mFBP and LSQI in-focus images comes at a price; namely, this fiber persist strongly to out-of-focus planes.

The relative intensity metric provides a means to interpret the fiber conspicuity in the DBT scanner images of the tilted ACR-DM phantom. For the present purpose, the metric is written $i(z,\phi)$ to allow for $z$-dependence and the rod length is set to $L=10.0$mm to match the ACR phantom fiber length. The in-focus plane images in Fig. 8 correspond to the relative intensity plots shown in Fig. 10. Compared to the shorter $1.9\text{\,}\mathrm{mm}$ rods of Fig. 5 the dependences in the relative intensity for the $10.0\text{\,}\mathrm{mm}$ rods is more exaggerated. In particular, the rod signals nearly parallel to the DXST have a much lower relative intensity than their shorter counterparts. This agrees with the visual assessment of the fiber signals in Fig. 8, where the 5^∘ fiber is difficult to see except when strong regularization is employed. The images in Fig. 9 shown for planes below the focus plane correspond to the curves plotted in Fig. 11. The relative intensity curves track the persistence of the small-angle ACR-DM fibers from Fig. 9 seen for regularized mFBP and LSQI. The corresponding relative intensity curves in Fig. 11 increase with greater regularization particularly for the low-angle rods reconstructed with mFBP and LSQI. This result demonstrates the trade-off between in-focus-plane orientation dependence of fiber conspicuity and depth blur of fibers nearly parallel to the DXST.

Up until the present section, we have focused on various aspects of the imaging of fiber-like structures embedded in uniform backgrounds. To investigate visualization of fiber-like signals in the presence of realistic background structure, we show an ROI from a clinical DBT case reconstructed by FBP, mFBP, LSQI, and LSQD. The chosen ROI has fibrous structures at different orientations so that some of the discussed trends in rod signal conspicuity can be illustrated.

In Fig. 12 an ROI from a clinical case is displayed in an arrangement similar to the ACR-DM images shown in Fig. 8; the same ROI is shown for FBP, mFBP, LSQI, and LSQD at the same regularization strengths used for the ACR-DM images. This particular ROI is chosen because of the presence of multiple fiber-like structures at various in-plane orientations; and the conspicuity of these structures can be compared among the algorithms and regularization strengths. In particular, we highlight a group of three nearly horizontal vessels that change in conspicuity dramatically depending on algorithm and regularization parameter setting. Their conspicuity varies in a manner consistent with the presented rod-signal simulations and ACR-DM phantom results.

The background texture in the patient images of Fig. 12 is seen to vary among the algorithms and parameter settings, and the differences are more obvious than what is seen in the ACR-DM images of Fig. 8. One aspect of the texture can be traced to imaging of fiber-like structures; the mFBP images show horizontal wisp artifacts for $\delta\geq 0.05$ and to a slightly lesser degree the same wispiness is seen in LSQI for $\gamma\geq-0.75$. These artifacts can be traced to vessels that are present at slices of different depths. For example, the indicated horizontal wisp in the mFBP image for $\delta=0.12$ can be traced to a vessel oriented parallel to the DXST in a plane that is 1.3 cm below the shown image. The horizontal wispiness is a direct result of the persistence of horizontal fiber-like structures from other slices at different depths.