1.

Introduction

Iterative image reconstruction (IIR) for x-ray computed tomography (CT) has been heavily investigated in recent years.¹ It is well known that IIR involves a substantially greater computational effort than direct reconstruction algorithms such as filtered back-projection (FBP). For CT, this issue is particularly acute because it is a high-resolution imaging modality, requiring the use of small voxels to achieve the inherent resolution of the CT data.² Because IIR is based on implicit estimation of the image, a complete representation of the subject within the scanner view domain is also necessary. Objects appearing in the x-ray projections such as the patient couch, which have no relevance to the imaging task, must also be accounted for in the reconstruction volume. Such a complete representation is only necessary for implicit methods and not direct reconstruction algorithms such as FBP, where the inversion equation need only be evaluated at points within an ROI. The combination of IIR requiring the representation of large volumes and use of small voxel exacerbates the computational burden issue for CT image reconstruction. To address this problem, we propose a method for direct ROI IIR, which will allow for the use of smaller reconstruction volumes focused on ROIs relevant to a given imaging task.

The main tool employed to realize ROI IIR appears in our conference record (CR) publication for the 2014 SPIE medical imaging conference.³ In the CR, we presented an offshoot of our adaptive steepest descent–projection onto convex sets algorithm for edge enhanced image reconstruction applied to clinical digital breast tomosynthesis (DBT) projection data. To reach the desired DBT image quality goals, we proposed to compare the estimated and DBT projection data after convolving with a derivative kernel. The resulting IIR algorithm was seen to behave analogously to a class of analytic image reconstruction methods called $\mathrm{\Lambda }$ or local tomography.^4–7 The resulting DBT images accentuated small structures relative to IIR employing the standard data fidelity term. Furthermore, the proposed derivative data fidelity yielded an IIR algorithm robust against heavily truncated projection data allowing for direct image reconstruction of a DBT ROI. It is this latter feature of the modified data fidelity term that we exploit for ROI IIR in CT.

The algorithm developed in Ref. 3 was designed to obtain a useful image after only a few iterations, but it was not run to convergence of the optimization problem upon which the algorithm was based. As a result, the obtained images depend on parameters characterizing the derivative filter in addition to other optimization and control parameters of the algorithm. In this article, we aim to characterize the solution of the proposed optimization problem and accordingly we employ an accurate solver developed by Chambolle and Pock (CP).^8,9 In this way, we narrow down the IIR parameter space to only those needed in specifying the optimization problem. With the understanding gained by thoroughly investigating the solution of the proposed optimization problem, design of IIR with severely truncated iteration is facilitated.

In Sec. 2, we present the proposed optimization problem and the CP algorithm derived to solve it. In Sec. 3, we characterize the solution of the optimization problem with ideal and noisy simulated data for sparse-view and standard CT acquisition. In Sec. 4, we focus on ROI image reconstruction, comparing the use of the proposed derivative weighted data fidelity term with that of a standard Euclidean data fidelity with no additional weighting. Finally, in Sec. 5, we demonstrate ROI reconstruction on actual CT scanner data.

2.

Optimization for CT IIR Using a Derivative Weighted Data Fidelity Term

We introduce an optimization problem designed to allow for direct ROI image reconstruction in CT, and present the algorithm used to solve the optimization problem. The CT data model is presented in Sec. 2.1; the particular data fidelity term of interest is introduced in Sec. 2.2, where the link with $\mathrm{\Lambda }$-tomography is also explained; incorporation of this fidelity into an optimization for CT IIR is discussed in Sec. 2.3; and finally, a pseudocode for the solver is given in Sec. 2.4.

2.3.

Optimization Problem for CT IIR

The optimization equation we investigate here is

Eq. (11)

$${\mathbf{f}}^{*}=\underset{\mathbf{f}}{\arg \min }\frac{1}{2}{\Vert (D_u+cI)(X\mathbf{f}-\mathbf{g})\Vert }_2^2\text{such that}\Vert (|\nabla \mathbf{f}|){\Vert }_1\le \gamma .$$

The operator $\nabla$ is a matrix representing a finite differencing approximation to the image gradient. The number of pixels/voxels in $\mathbf{f}$ is taken to be $N$ and the spatial gradient operator $\nabla$ yields a vector-valued image, e.g., $\overrightarrow{\mathbf{z}}$, of size $2N$ or $3N$ depending on whether the image is 2-D or 3-D, respectively. The vector symbol is used to indicate a spatial vector-valued image, and the spatial magnitude operator, $|·|$, converts a spatial vector-valued image to a scalar image by taking the vector magnitude at each pixel/voxel. The quantity $\nabla \mathbf{f}$ is the spatial gradient of $\mathbf{f}$; $|\nabla \mathbf{f}|$ is the gradient magnitude image (GMI); and ${\Vert (|\nabla \mathbf{f}|)\Vert }_1$, the sum over the GMI, is the image total variation (TV). The optimization equation, Eq. (11), seeks the image with lowest combined derivative and Euclidean data discrepancy amongst all images with TV less than or equal to $\gamma$.

This particular optimization problem is one variation of many that combines some form of data fidelity with a term involving image TV. The TV term is useful for exploiting sparsity in the GMI for image reconstruction problems with limited data. We use it exactly for this reason; the direct ROI reconstruction problem is a form of limited data acquisition, because this problem is equivalent to the projection truncation problem for IIR. The reason for this is that in the implicit solution of Eq. (2), only transmission rays that pass through the ROI are used. This contrasts with direct reconstruction through FBP, where all projection data are used in filtering before back-projecting to the ROI. The particular form of the gradient sparsity exploiting optimization problem, where the constraint is placed in the image TV, is particularly convenient because the data fidelity metric is actually the sum of two fidelity metrics. Furthermore, for computer simulation studies, as performed in this article, the phantom TV provides a good reference for selecting the parameter $\gamma$.

To completely specify Eq. (11) for the following CT simulations, we select a particular form of $D_u$. As required, $D_u$ is chosen to be purely antisymmetric and it represents convolution with a 21-point kernel. This kernel is generated by smoothing an antisymmetric differencing kernel with a Gaussian of standard deviation $\omega$ measured in detector bin units. The smoothing of the finite differencing kernel is useful for controlling noise texture in the reconstructed images. The kernel is shown in Fig. 1 for different values of $\omega$.

Fig. 1

Kernel used in specifying the derivative filter $D_u$ for different values of $\omega$, the standard deviation of the Gaussian smoothing.

Aside from the usual scanning geometry parameters, there are three parameters—$\gamma$, $c$, and $\omega$—involved in the optimization problem of study in Eq. (11). The goal of this article is to demonstrate the effect of each of these parameters on the solution of Eq. (11) in order to gain intuition about this optimization problem for various CT imaging configurations and in particular for direct ROI IIR. For real data applications, such as the DBT reconstructions presented in Ref. 3, an IIR algorithm with a truncated iteration may be used. In this case, we do not expect to achieve the solution to Eq. (11) and as a result, the attained reconstructed volume depends on all the parameters of the IIR algorithm further complicating the characterization of the algorithm. Thus, before a thorough investigation of truncated IIR motivated by Eq. (11), it is important to characterize the solution of this optimization problem. To this end, we employ the CP algorithm for its accurate solution.

3.

Characterization with Simulated Untruncated Projection Data

The simulated 2-D CT data are based on the breast CT imaging application. The computer generated phantom shown in Fig. 2 models fat and fibro-glandular tissues. The complex structure is obtained by generating an image with power-law filtered white noise, followed by thresholding. The structure model is described in Ref. 13. The scanning geometry is 2-D fan-beam CT with full circular coverage, as what would be available in the midplane of an actual breast CT scan. The source-to-isocenter distance is 36 cm, and the source-to-detector distance is 72 cm. A linear detector with 1024 bins is modeled, and for the majority of the simulations, 256 projection views are taken. The exceptions are 1024 views are used for the studies without TV regularization, and 64 views are employed in one study examining the possibility of sparse-view sampling. In this initial series of characterization studies, we examine image reconstruction from ideal noiseless data, where the projections are generated from a phantom defined on the same $512\times 512\text{pixel}$ grid used for the image reconstruction. We also perform studies for inconsistent CT data with noise modeling a typical breast CT scan, and projection operating on mismatched grids. The phantom projections are generated from a $1024\times 1024\text{pixel}$ grid and the image estimate projections are generated from a $512\times 512\text{pixel}$ grid.

Fig. 2

Breast phantom for computed tomography (CT) and its corresponding gradient magnitude image (GMI). (a) The linear attenuation map of the phantom in the gray scale window $[0.174,0.253]{\mathrm{cm}}^{-1}$. The attenuation values for the tissues, fat, and fibro-glandular, simulated in the phantom are 0.194 and $0.233{\mathrm{cm}}^{-1}$, respectively. These values model a monochromatic x-ray source at 50 keV. (b) The GMI $[0.0,0.1]{\mathrm{cm}}^{-1}$ illustrating the test phantom’s sparsity in the GMI.

Comparison of images presented in this work may not be straight-forward because of the loss of absolute gray level information. For some series of images, in order to have objective image comparisons, we compute the mean, $f_{\mathrm{mean}}$, and standard deviation, $f_{\mathrm{std}}$, of the pixel values within a specified region contained just inside the boundaries of the test phantom or the field-of-view (FOV), whichever is smaller. The gray scale range is computed by

To assist in the interpretation of the CT IIR results, we show FBP and $\mathrm{\Lambda }$-tomography images in Fig. 3. The $\mathrm{\Lambda }$-tomography image shows the expected enhancement at tissue borders and loss of the absolute gray level. Again, the reason $\mathrm{\Lambda }$-tomography might be an interesting choice in addition to the edge enhancement is that the filter is local, allowing for projection truncation.

Fig. 3

Analytic reconstruction from noiseless data comprised of 1024 projections. (a) Reconstruction by filtered back-projection (FBP), using a pure ramp filter, and (b) $\mathrm{\Lambda }$-tomography where the data are convolved with a finite differencing approximation to the second derivative prior to back-projection. The $\mathrm{\Lambda }$-tomography image aids in the interpretation of the images at low iteration numbers of the presented iterative image reconstruction (IIR) algorithm. The gray scale is computed by use of Eq. (15).

The first set of IIR images result from a CP algorithm designed to minimize the standard quadratic data fidelity ${\varphi }_{\mathrm{Eucl}}(\mathbf{f})$. The pseudocode appears in Algorithm 2 of Ref. 9. Intermediate images of this algorithm are shown in Fig. 4, where the 1024-view data are noiseless and perfectly consistent. The behavior of the progression of images, starting from a zero initial image, is generic to other first-order nonsequential algorithms applied to minimization of ${\varphi }_{\mathrm{Eucl}}(\mathbf{f})$. The image after the first iteration is blurry; in fact, by a similar argument to the one presented in Sec. 2.2, one can show that this image is approximately the result of back-projection applied to the projection data. As the iterations progress, the images become less blurry, and eventually the test phantom is obtained.

Fig. 4

Iterative image reconstruction from noiseless projection data containing 1024 views by use of a Chambolle and Pock (CP) algorithm designed to minimize ${\varphi }_{\mathrm{Eucl}}(\mathbf{f})$. The image sequence shows intermediate iterations, and the iteration number is indicated in each panel. The gray scale window appears above each image. The gray scale is computed by use of Eq. (15).

The next set of IIR images result from a CP algorithm designed to minimize the proposed data fidelity ${\varphi }_{\mathrm{combo}}(\mathbf{f})$ with $\omega =0$ and $c=0$. The same pseudocode from Algorithm 2 of Ref. 9 can be employed to minimize this data fidelity. Similar to gradient descent, one can show that the first iteration of this CP algorithm is proportional to a $\mathrm{\Lambda }$-tomography image if a zero initial estimate is used. Intermediate images of this algorithm are shown in Fig. 5, where the 1024-view data are noiseless and perfectly consistent. Indeed, we observe that the first panel shows an image proportional to that of $\mathrm{\Lambda }$-tomography and as the algorithm proceeds to convergence we observe that the images tend toward the phantom. Minimization of both ${\varphi }_{\mathrm{Eucl}}(\mathbf{f})$ and ${\varphi }_{\mathrm{combo}}(\mathbf{f})$ results in accurate recovery of the test phantom, but the approach to this solution is quite different; the first iterate in minimizing ${\varphi }_{\mathrm{combo}}(\mathbf{f})$ already has much of the tissue structure information relative to that of minimizing ${\varphi }_{\mathrm{Eucl}}(\mathbf{f})$. Although the image appearance of the early iterations in Fig. 5 has been explained, it is not obvious that the phantom should be recovered accurately at convergence—in light of the finite-differencing filter applied to both system matrix and data in ${\varphi }_{\mathrm{combo}}(\mathbf{f})$.

Fig. 5

Iterative image reconstruction from noiseless projection data containing 1024 views by use of a CP algorithm designed to minimize ${\varphi }_{\mathrm{combo}}(\mathbf{f})$. The image sequence shows intermediate iterations, and the iteration number is indicated in each panel. The gray scale window appears above each image. The gray scale is computed by use of Eq. (15).

In order to understand the phantom recovery, we apply singular value decomposition (SVD) in a manner similar to that of Ref. 14. We note that under ideal noiseless conditions, the minimizer of ${\varphi }_{\mathrm{combo}}(\mathbf{f})$ with $\omega =0$ and $c=0$ satisfies the following linear equation:

Fig. 6

Maximum and minimum singular values for the combined system matrix $D_{u}X$ for a two-dimensional CT configuration with image pixel array $N\times N$ and projection data set consisting of $2N$ projections over a $2\pi$ scanning arc and a linear detector array with $2N$ detection elements. The solid circles represent singular values arrived at through direct SVD of a small system with $N$ ranging from 18 to 64. The solid diamond represents extrapolation to the case of interest at $N=512$. The linear fitting model is displayed in the legend of each graph.

It is also interesting to observe the evolution of the CP algorithm iterates when the sampling ratio is more challenging. Reducing the view number to 256 projections results in a system matrix $D_{u}X$ with a much larger condition number. Intermediate images of the CP algorithm are shown in Fig. 7, where the 256 view data are noiseless and perfectly consistent. We again observe that the first panel shows an image proportional to that of $\mathrm{\Lambda }$-tomography and as the algorithm proceeds to convergence we observe that the images seem to tend toward the phantom; however, the converged image shows noise-like texture in the images resulting from the insufficiency of the number of projections. Interestingly, however, the absolute gray level appears to be accurately recovered despite the poor overall conditioning of $D_{u}X$ with 256 projections.

Fig. 7

Same as Fig. 5 except that the projection data contain only 256 views. The gray scale is computed by use of Eq. (15).

As discussed in Sec. 2.3, the view angle undersampling can be dealt with by exploiting GMI sparsity. To demonstrate this, we solve the optimization problem in Eq. (12) using the pseudocode in Algorithm 1 and again model noiseless consistent 256-view projection data. Because the data are ideal, the parameters $c$ and $\omega$ are again set to zero. The TV constraint parameter is selected to be that of the phantom, $\gamma ={\gamma }_0$. The resulting sequence of images is shown in Fig. 8, and they resemble those of Fig. 7 except that at iteration 10 and after the noise-like artifacts seem to be removed by the addition steps constraining the image TV. Furthermore, at convergence, the phantom is recovered to high accuracy. The accurate phantom recovery for this example, of course, depends on the fact that we know the correct TV constraint value $\gamma$, but the purpose here is to show that an ideal program can recover the phantom under ideal conditions. That it can do so is not obvious because the data discrepancy objective is different than the standard Euclidean form. Although the $\mathrm{\Lambda }$-tomography interpretation of the low iteration images is useful, it is important to keep in mind that the specific evolution of the images depends on the values of all of the parameters in Algorithm 1 in addition to parameters of the optimization problem Eq. (12). The converged image, on the other hand, depends only on parameters of the optimization problem.

Fig. 8

Same as Fig. 7 except that the image sequence results from use of Algorithm 1, which is designed to solve total variation (TV)-constrained minimization of ${\varphi }_{\mathrm{combo}}(\mathbf{f})$, see Eq. (12). Note that the first images resemble that of $\mathrm{\Lambda }$-tomography. We also observe that at convergence the phantom is recovered accurately. The gray scale is computed by use of Eq. (15).

In order to further test the robustness against view angle undersampling, we reduce the ideal data set to only 64 projections. The resulting image sequence appears in Fig. 9. We observe that Algorithm 1 can accurately recover the test phantom as the image estimates near convergence. The $\mathrm{\Lambda }$-tomography interpretation of the low iteration images, however, suffers due to the streaks from the view angle undersampling. Nevertheless, the insensitivity of the proposed data discrepancy term to absolute gray level does not appear to hinder the ability of Algorithm 1 to recover the phantom from sparse view data.

Fig. 9

Same as Fig. 8 except that the simulated projection data contain only 64 views. In this case, the low iteration images are contaminated by view-angle undersampling artifacts, which appear in the form of streaks. As the image estimates approach convergence, we observe that the phantom is again accurately recovered. The gray scale is computed by use of Eq. (15).

Finally, we investigate the response of Algorithm 1 to inconsistent projection data. Two forms of inconsistency are introduced in the data model. First, there is mismatch between the discrete phantom grid, $1024\times 1024$, and reconstructed image grid, $512\times 512$. Second, noise is introduced with a Poisson-like Gaussian distribution, i.e., the variance of the Gaussian is set equal to the mean, where the mean is the average transmission at each detector bin assuming $7.5\times {10}^4$ photons are incident to each bin prior to passing through the subject. Thus, we are modeling a low-dose CT scan as what is typically proposed for breast CT. Images corresponding to the accurate solution of Eq. (12) with the TV constraint parameter selected to be that of the phantom, $\gamma ={\gamma }_0$, and various values of the parameters $\omega$ and $c$, are displayed in Fig. 10. The parameter $\omega$ clearly impacts noise texture, and here, we only aim to show how image quality changes with this parameter. Optimal setting of $\omega$ can only be determined based on the particular image task of the CT scan. Although $c$ has little impact on the reconstructed images for these simulations with untruncated projections and full image representation, this parameter helps to control gray level for ROI imaging.

Fig. 10

Images reconstructed from 256 projections with an inconsistent data model by Algorithm 1. In the top row, the derivative smoothing width is fixed at $\omega =1$, and for the bottom row the combination parameter is set to $c=0$. The gray scale for each of the images is [0.171, 0.233]. As can be seen $c$ has no perceptible effect on these images where the data contain no truncated projections. The derivative smoothing width, on the other hand, has a clear impact on the noise texture.

4.

Characterization for CT ROI Imaging

In this section, the CT simulations aim at illustrating the utility of Algorithm 1 for ROI imaging, where an incomplete image representation is used that spans only a given ROI. This study is performed in two steps: first, using a full image representation, we reconstruct the FOV of a truncated projection acquisition; and second, using only image pixels in this same FOV, we perform image reconstruction. By the use of the term FOV, we specifically refer to the collection of all pixels that are illuminated by all views in the data set. We use this term instead of “interior tomography,”¹⁵ which has a specific meaning in 2-D image reconstruction where the aim is to recover a continuous region interior to the support of the subject from projections illuminating only this interior region. Furthermore, for interior tomography, the data model is the continuous 2-D Radon or x-ray transform, while for the present FOV imaging, the data model is the discrete x-ray transform. We make the distinction between ROI and FOV in order to allow for the design of a scan where the FOV is larger than the ROI in order to avoid potential artifacts at the FOV border.

First, we investigate the impact of truncated data while keeping the image representation complete, i.e., the image pixels cover the region where the object support is nonzero. The simulated projection data are ideal and noiseless, containing 256 projections. Each of the projections is truncated so that the FOV of the scan is a centered circle of diameter 9 cm, i.e., half the width of the image array. The left column of images in Fig. 11 corresponds to converged solutions of Eq. (12) for $\omega =0$ and $c=0$, and for comparison, the right column of images results from the same optimization problem without the proposed weighting $D_u$. For this configuration, it is not clear what is the optimal choice of the TV constraint parameter $\gamma$, even though we know what is the test phantom. Using the TV of the phantom within the FOV, ${\gamma }_{\mathrm{FOV}}$ and the TV of the entire phantom ${\gamma }_0$ as lower and upper bounds on $\gamma$, respectively, we reconstruct images for four intermediate values of $\gamma$ parametrized by ${\gamma }_w$

Fig. 11

Images reconstructed from truncated, ideal noiseless data. The shown images on the left column are converged solutions to Eq. (12) with $c=0$ and $\omega =0$, and those on the right column are converged solutions of the same optimization without the weighting $D_u$, i.e., TV-constrained Euclidean data discrepancy minimization. The value of ${\gamma }_w$ indicated in each panel determines the TV constraint parameter $\gamma$ by Eq. (17). The horizontal lines indicate the location of the image profiles plotted in Fig. 12. The gray scale is set to [0.165, 0.265].

For ${\gamma }_w=1$, use of the weighting $D_u$ seems to make little difference; both images show cupping although that of the $D_u$ weighting image has slightly less cupping than when no weighting is used. As ${\gamma }_w$ decreases, both sets of images show a decrease in resolution particularly in going from ${\gamma }_w=0.5$ to 0.25; however, for the $D_u$ weighting images edges within the FOV appear to remain sharp relative to the edges in the images generated without the proposed weighting. Moreover, the $D_u$ weighting images show markedly less cupping. A quantitative sense of these observations is obtained in the profile plots of Fig. 12.

Fig. 12

Profiles of the images shown in Fig. 11 compared with the phantom profile.

Second, we investigate the impact of employing an incomplete image representation, where pixels cover only the region of the FOV from the previous study, a centered 9-cm diameter circular ROI. Note that in this case it does not matter whether the projections are truncated to this 9 cm ROI or not. For untruncated projections, measurements corresponding to the rays not intersecting the ROI end up not playing any role in the IIR algorithm. Converged images—with and without the $D_u$ weighting—obtained by Algorithm 1 appear in Fig. 13. It is in the use of the smaller ROI representation where we see a possible major advantage of the $D_u$ weighting. The $D_u$ weighting reduces the data discrepancy incurred by using too small an image representation. As a result, the TV constraint can be applied without destroying structure information in the ROI while the use of the TV constraint without the $D_u$ weighting destroys this information. Inspecting the gray scale windows for both plots, we note that the automatically computed gray scale for the $D_u$ weighting image is substantially lower than the values of the actual phantom, while the image without weighting has approximately the correct absolute gray scale (aside from the fact that the structure information is lost). That the absolute gray scale of the $D_u$ weighting is lost is the primary motivation for introducing the parameter $c$.

Fig. 13

Images reconstructed from ideal noiseless data containing 256 projections. The image representation is restricted so that all pixels are contained in a 9-cm diameter ROI. The left image is a converged solution to Eq. (12) with $c=0$ and $\omega =0$, and the the right image is a converged solution of the same optimization without the weighting $D_u$, i.e., TV-constrained Euclidean data discrepancy minimization. The TV constraint parameter is set to the TV of the phantom within the ROI. The gray scale is computed by use of Eq. (15).

For the final simulations, we again employ the ROI-only representation, but model inconsistent projection data with grid mismatch and measurement noise, corresponding to $7.5\times {10}^4$ photons incident to each detector bin. Again, the data acquisition configuration is chosen so that the FOV coincides with the 9-cm diameter ROI. In the grid of images shown in Figs. 14 and 15, results are shown for various values of $\gamma$ and $c$, while $\omega$ is fixed to one. The value of $\gamma$ is parametrized in these figures with the multiplicative factor ${\gamma }_f$

Fig. 14

Images reconstructed from inconsistent noisy data containing 256 projections. Each of the projections is truncated so that the field-of-view (FOV) of the scan is a centered circle of diameter 9 cm and the image representation is restricted so that all pixels are contained in the FOV. The array of images results from accurate solution to Eq. (12) for $\omega =1$ and various values of $c$ and the TV constraint parameter is set according to Eq. (18). The values of ${\gamma }_f$, see Eq. (18), indicated in the left most images remain constant for each row, and the values of $c$ indicated in the bottom row of images remain constant for each column. The image gray scale is [0, 0.35], wide enough to include all images for the sake of comparison. For individual images, of course, a narrower gray scale may be employed to improve contrast. The test phantom is displayed in the lower, right panel for reference.

Fig. 15

Same as Fig. 14 except that the gray scales are determined by Eq. (15) and written above each of the corresponding image panels.

The leftmost column of images in Fig. 15 shows results for $c=0$. Each of these images show the structural information of the phantom within the FOV, but the gray values are quite low and independent of $\gamma$. Within this column, the TV constraint is an effective control over image regularization, and the image corresponding to ${\gamma }_f=1$ shows almost no noise artifacts while suffering minimal blurring due to the regularization. As $c$ increases, the gray levels of the reconstructed images also increase. The restoration of the gray levels, however, comes at the price of the TV constraint becoming less effective. Following the ${\gamma }_f=1$ row, the increasing data discrepancy caused by increasing $c$ results in image details being smoothed away. Inspecting the $c=0.05$ column, the image details can be recovered by loosening the TV constraint, but there is a clear loss of detail relative to the $c=0$ images. Nevertheless, for some imaging applications, it may be important to employ the $c$ parameter to recover absolute gray level information.

5.

Application of ROI IIR to Severely Truncated Projections from Actual CT Scanner Data

To demonstrate application of the proposed algorithm to actual scanner data, we employ an XCounter CT scan of a rat. The data set contains 1000 projections with a $768\times 512$ detector array where the detection elements are 0.1 mm in width. In order to test the ROI capability, we truncate the projection data to include only transmission rays passing through a cylindrical subvolume of radius 8 mm and height 18 mm. This volume contains part of the rat jaw and a fore paw. For reference, we employed FBP to reconstruct the complete untruncated data set, and display slices that traverse the chosen ROI in Fig. 16.

Fig. 16

From left to right are a transaxial slice, normal to the $z$-axis, and vertical slices normal to the $y$- and $x$-axes reconstructed from the XCounter CT scan of a rat by a standard FBP implementation. The display plane resolution is 0.1 mm. In order to reduce noise, three 0.1-mm thick slices along the view axis are averaged. The display gray scale is $[-0.08,0.45]$. The boxed regions indicate the display bounds of Fig. 17.

Returning to the truncated data set, we apply FBP, $\mathrm{\Lambda }$-tomography, and the proposed ROI IIR algorithm with $c=\omega =0$ for two different values of the TV constraint, $\gamma$. The constraint values correspond to an average gradient magnitude of $1.25\times {10}^{-3}({\mathrm{cm}}^{-1})$ and $6.0\times {10}^{-4}({\mathrm{cm}}^{-1})$. The resulting images are shown in Fig. 17. The FBP images show the usual cupping in the FOV of the scan, and it is difficult to see structure near or beyond the FOV. The images generated by $\mathrm{\Lambda }$-tomography do not show a discontinuity at the FOV border, and the edges of structures are clearly seen but the gray levels are, as expected, lost. The ROI IIR algorithm appears to retain some gray level information without suffering from strong discontinuities at the FOV border. The difference in behavior near the FOV border is most striking in the images of the middle column of Fig. 17 which show the rat fore paw. We point out that this study on CT scanner data is only a preliminary demonstration and we did not attempt to optimize image quality over the parameters $c$, $\omega$, and $\gamma$. Likewise, the implementations of FBP and $\mathrm{\Lambda }$-tomography are basic because these images serve only as a reference. For the FBP images, we did not explore alternate filter designs nor attempt to counteract the cupping artifact by postprocessing or extrapolating projections.

Fig. 17

ROIs reconstructed by, from bottom row to top row, FBP, $\mathrm{\Lambda }$-tomography, ROI IIR with $c=\omega =0$ and $\gamma ={\gamma }_1$ and same for $\gamma ={\gamma }_2$. The constraint parameters ${\gamma }_1$ and ${\gamma }_2$ correspond to average gradient magnitudes of $1.25\times {10}^{-3}({\mathrm{cm}}^{-1})$ and $6.0\times {10}^{-4}({\mathrm{cm}}^{-1})$, respectively. From left to right are a transaxial slice, normal to the $z$-axis, and vertical slices normal to the $y$- and $x$-axes. The display plane resolution is 0.1 mm. In order to reduce noise, three 0.1-mm thick slices along the view axis are averaged. The display gray scale that was set by eye is given above each image.

6.

Discussion

We have presented an IIR algorithm based on an alternative data fidelity term that is designed to facilitate ROI reconstruction. The use of the restricted image representation can allow for IIR with high-resolution image representation without incurring prohibitively high computational and memory demands. Furthermore, the use of the proposed filtering on the data discrepancy preserves edges allowing for a high degree of regularization before structures in the subject are blurred away.

The present investigation mainly focuses on image properties of the converged solution of TV constrained data discrepancy minimization, where the images depend on three parameters: $\gamma$ controlling regularization, $\omega$ affecting noise texture, and $c$ allowing for gray level recovery. It may be that for practical application of the proposed data fidelity term, IIR algorithms operate with truncated iteration. Characterizing such algorithms becomes more complicated as other algorithm parameters in addition to the three of the optimization problems affect the reconstructed image when convergence of the corresponding optimization problem is not achieved. The established link with $\mathrm{\Lambda }$-tomography, however, can aid in the design of IIR using the proposed derivative weight and operating with truncated iteration.