Real-time out-of-plane artifact subtraction tomosynthesis imaging using prior CT for scanning beam digital x-ray system

Purpose: The scanning beam digital x-ray system (SBDX) is an inverse geometry fluoroscopic system with high dose e ffi ciency and the ability to perform continuous real-time tomosynthesis in multiple planes. This system could be used for image guidance during lung nodule biopsy. However, the reconstructed images su ff er from strong out-of-plane artifact due to the small tomographic angle of the system. Methods: The authors propose an out-of-plane artifact subtraction tomosynthesis (OPAST) algorithm that utilizes a prior CT volume to augment the run-time image processing. A blur-and-add (BAA) analytical model, derived from the project-to-backproject physical model, permits the generation of tomosynthesis images that are a good approximation to the shift-and-add (SAA) reconstructed image. A computationally practical algorithm is proposed to simulate images and out-of-plane artifacts from patient-specific prior CT volumes using the BAA model. A 3D image registration algorithm to align the simulated and reconstructed images is described. The accuracy of the BAA analytical model and the OPAST algorithm was evaluated using three lung cancer patients’ CT data. The OPAST and image registration algorithms were also tested with added nonrigid respiratory motions. Results: Image similarity measurements, including the correlation coe ffi cient, mean squared error, and structural similarity index, indicated that the BAA model is very accurate in simulating the SAA images from the prior CT for the SBDX system. The shift-variant e ff ect of the BAA model can be ignored when the shifts between SBDX images and CT volumes are within ± 10 mm in the x and y directions. The nodule visibility and depth resolution are improved by subtracting simulated artifacts from the reconstructions. The image registration and OPAST are robust in the presence of added respiratory motions. The dominant artifacts in the subtraction images are caused by the mismatches between the real object and the prior CT volume. Conclusions: Their proposed prior CT-augmented OPAST reconstruction algorithm improves lung nodule visibility and depth resolution for the SBDX system. C 2014 American Association of Physicists in Medicine. [http: // dx.doi.org / 10.1118 / 1.4896818]


INTRODUCTION
Fluoroscopy is poor at detecting lung nodules that are ≤20 mm in diameter and may miss tumors that are concealed by hilar vessels, mediastina, the heart, or the diaphragm. 1,2 Digital tomosynthesis is an imaging technique that uses radiation dose levels that are lower than those of standard CT fluoroscopy while at the same time providing depth information that is not available in radiographic fluoroscopy. 3,4 Digital tomosynthesis uses a series of projections from multiple angles to generate slice images with high in-plane (x-y) resolution. We aim to use tomosynthesis to replace fluoroscopy as a real-time image guidance tool for lung nodule transbronchial needle biopsy (TBNbx).
A real-time tomosynthesis approach is available based on the scanning beam digital x-ray (SBDX) hardware, originally developed for cardiac applications by NovaRay, Inc. 5 The SBDX system is an inverse geometry system with high dose efficiency and the ability to perform continuous realtime tomosynthesis. 6 The scanning source obviates the need for mechanical motion, permitting acquisition of digital tomosynthesis data sets at speeds limited only by the detector readout rate. In the SBDX system, the angle of the xrays passing through a voxel (tomographic angle) is approximately 5 • -15 • . The SBDX system has been investigated for lung cancer detection, 3D tracking of high-contrast objects in cardiac catheter tracking, and measurement of vessel dimensions at various magnifications and angulations without calibration. [7][8][9] Our application for the SBDX system aims to provide real-time and accurate visualization of the biopsy needle and the lung tumor nodule at the moment for biopsy. This image-guidance goal is not met by the current image reconstruction algorithms.
Shift-and-add (SAA) is a well-known analytical image reconstruction algorithm for digital tomosynthesis, 10 which is also used on the SBDX system. 6 When parallel motion of the tube and/or detector is implemented, the SAA algorithm involves shifting each of the projection images by a specific amount and then adding them together. By selecting the shifts-per-plane amount correctly, objects in a given plane can be brought into sharp focus. In addition to SAA, two deblurring algorithms that have received significant attention in recent years are matrix inversion tomosynthesis (MITS) and filtered backprojection (FBP). 11,12 However, the relative motion between the x-ray source and the detector covers a limited angular range; therefore, the depth (z) resolution is much lower than the in-plane resolution. The structures outside the focused plane are blurred and superimposed onto the reconstructed slices. The slice sensitivity profile for a weighted filtered backprojection tomosynthesis algorithm is inversely proportional to the tangent of the half scan angle. 13 Therefore, additional information (e.g., prior CT) is needed to augment the reconstruction to improve the depth resolution of the digital tomosynthesis.
In most TBNbx procedures, 3D thorax CT images are used to plan the biopsy paths. 14 The patient-specific prior CT volume can also be used to augment the tomosynthesis image reconstructions and processing to suppress the out-of-plane structures. Because the TBNbx procedure is performed with a prior CT, the CT volume (and its projections) is readily available. We want to utilize the patient-specific prior CT volume to augment the run-time reconstruction of the SBDX system. Several previous studies have used available prior knowledge to improve the depth resolution in digital tomosynthesis images. These approaches fall into three categories: registrationbased, data-filling-based, and subtraction-based approaches.
The registration-based approaches estimate motions and deformations of the reference CT volume using projectionsto-volume registration. [15][16][17] A deformation map is calculated that minimizes the differences between the synthesized projection images and the detected projections. Recent work focuses on using known motion models, free-form deformation, optical flow, etc., with limited-angle or sparse projections to reduce the motion artifact in cone-beam CT and tomosynthesis. [18][19][20] The advantages of registration-based approaches include well-developed nonrigid transformation models and relatively mature volume-to-projection registration solvers. Thus, registration-based algorithms can have very high accuracy even for systems with small tomographic angles. However, most of the registration-based algorithms require iterative solvers for the registration, which are usually slow and limit their use in real-time systems. The registration algorithms often require multiple projection/backprojection steps that have high computational demands.
The data-filling-based approaches use a low-resolution CT scan followed by a high-resolution tomosynthesis scan and then combines both projections for CT reconstructions. 21 The low-resolution CT volume augmentation is equivalent to filling in the (nonacquired) empty Fourier space of the tomosynthesis data with some smooth and meaningful information. The data-filling-based methods require highly accurate alignment between the CT and tomosynthesis scans. Any motions and deformations between the two systems may require proper compensation. In Zeng et al.'s work, an iterative algorithm is also used to combine the CT and tomosynthesis data. In addition, the low-resolution CT projections may not improve the depth resolution of the image details.
The subtraction-based approach proposed by Kim et al. subtracts projection components stemming from outside the region-of-interest (ROI) from the measured projection data. 22 Then, the preprocessed projections used for reconstruction only consider the structures within the ROI. The subtraction can suppress the artifacts caused by the out-of-ROI structures. However, the selected ROI in this approach is a 3D volume that includes multiple slices. Thus, the subtraction does not reduce the out-of-plane artifacts within the ROI. The depth resolution is not improved by the subtraction in the projection space. This approach also relies on accurate alignment in the projection space, which is a challenging and active area of research. 18,20,23 Yoon et al. proposed a fast image-based SBDX system analytical simulation model that produces equivalent tomosynthesis SAA reconstructions from a CT scan. 24,25 In the analytical blur-and-add (BAA) model, tomographic blurring of the SBDX geometry is modeled as a rectangular blur function, whose size varies as the distance from the focal plane times the tangent of the half-tomographic angle. Summation through the blurred slices provides an accurate approximation of the SAA tomosynthesis reconstruction. An important advantage of the BAA model is that each simulated slice is a summation of the blurred planes inside the entire volume. Therefore, the blurred out-of-plane structures can be easily computed by summing the selected blurred planes. Subtracting these simulated artifacts can be used to improve the depth resolution in the reconstruction. In this paper, we introduce a subtraction-based artifact suppression method, out-of-plane artifact subtraction tomosynthesis (OPAST), for the SBDX system. The BAA analytical model is modified to simulate images with an extended field-of-view (FOV) using the prior CT volume. After finding a good alignment between the SAA reconstruction and the BAA simulation, subtraction is used to suppress the out-of-plane artifacts. Figure 1 presents a schematic of the proposed OPAST algorithm.

2.A. SBDX system
The SBDX system consists of a scanning beam x-ray source and a 10.8 × 5.5 cm CdZnTe detector array, which are mounted on a C-arm gantry. 7 The scanning source consists of a 44-mm-thick, 23 × 23 cm collimator, which has an array of 100 × 100 holes with a pitch of 2.3 mm. The CdZnTe detector comprises 32, 13.5 × 7.5 mm tiles arranged in an 8 × 4 array. Each tile has 12 × 12 photon counting elements on a 1.14 mm pitch. The detector can perform at count rates in excess of 2 × 10 6 counts/s mm 2 with small cross talk between pixels. We changed the source-to-detector distance (D sd ) to 100 cm and the central-slice-to-source distance (D s0 ) to 60 cm for which the tomographic angles of the ROI are 15 • -20 • and F. 1. Schematic of the proposed OPAST algorithm for the SBDX system. The SAA reconstruction is first registered with the extended field-of-view BAA simulated image (indicated by the rectangles) and then subtracted using the BAA simulated out-of-plane artifacts. 5 • -15 • in the x and y directions, respectively. 26 The x-ray source can scan the 100 × 100 focal points 60 times in 1 s. For each reconstructed frame, the detector data are then the sum of several detector reads for each scanned focal spot. The number of photon counts per pixel typically averages 0.7 per capture period (1.28 µs).
The conventional SBDX image reconstruction directly applies the SAA algorithm to the detector photon count data followed by a pattern correction technique. 7 The detector data are usually zeros or small numbers, which is the main reason for using photon counts directly in the reconstruction instead of logarithmic data. The logarithmic data, commonly used in CT reconstructions, are also desired for the SBDX system in CT-augmented reconstruction. To obtain more photon counts per pixel for logarithmic computations, we propose to increase the detector pixel size by binning the detector to 48 × 24 with 2.28 × 2.28 mm 2 as well as by applying temporal averaging (i.e., reducing the frame rate). We also suggest downsampling the source plane to 50 × 50 points, such that a 4x increase in data-per-focal spot can be acquired at the same frame rate. Reconstructions are provided at lower frame rates of 5-10 fps, which results in total photon counts above 30 per pixel per reconstruction, and the photon statistics can then be approximated as a Poisson random variable.
The 2D system geometry and coordinates are shown in Fig. 2. The source focal spots lie on the line s parallel to the x axis. The flat panel detector pixels lie on the line u parallel to the x axis. The planes of interest in the volume are parallel to the source and detector planes. The centers of the source plane and detector plane align vertically along the z axis. The origin of z is defined at the source plane. D sd denotes the source-to-detector distance, and D s0 denotes the sourceto-central slice distance. For a source-to-reconstruction-plane distance z, the width of a back projected detector pixel is F. 2. Schematic 2D representation of the geometry and coordinates of the SBDX system.
where ∆ u denotes the detector pixel width. The shift of a voxel between two consecutive projections is where ∆ s denotes the pitch of the source focal spot.

2.B. SAA image reconstruction
The distance-driven (DD) method 27 is used as a discrete version of projection and backprojection for the SBDX system. The DD method uses a rectangular approximation of the footprints of a voxel on the detector. This method can be implemented with efficient memory access by separating and caching the footprint weights in the x and y directions. Mathematically, the DD projectors can be represented in the general form as where f (x,z) andf (x,z) denote the original image and the backprojected image at 2D spatial location x and z, respectively, and g(s,u) denotes the projection at pixel u from the source focal spot position s. The function F(s,u; x,z) denotes the value of the rectangular footprint of the voxel centered at (x,z) where ∆ z denotes the slice thickness, and x shift is The rectangular function in this paper is defined as The parameter l θ (s,u) is the amplitude function that approximates the ratio of the length of a ray passing through a slice to the slice thickness In our proposed 100 cm D sd and 60 cm D s0 SBDX system setup, we also use 2 × 2 detector pixel binning and 2 × 2 downsampling in source focal points to obtain more photon counts per pixels with the same frame rate. By summing more repeated x-ray passes together, the detector data (photon count) are high enough that the logarithm can be executed without creating large errors. Then, the SAA reconstructs an image of accumulated linear attenuation coefficients f recon (x,z) where y(s,u) denotes the logged detector data. The image voxels are reconstructed using many beamlets, and the number of overlapped beamlets per voxel is not uniform. The backprojection image needs to be normalized by the beamlet overlapping weights [the denominator in Eq. (8)].

2.C. BAA simulation model
Assuming the centers of the detector and source planes are perfectly aligned, lines connecting opposite edges of the detector and source planes intersect at where l det and l src denote the widths of the detector and source planes, respectively. We called the distance from the intersect point to the source plane the critical distance (D c ) because it divides the area between the source and detector into two regions. Most of the voxels above D c (i.e., point A in Fig. 2) can "see" the entire source but only part of detector. The x-ray beam, passing through this region, is mainly truncated by the finite size of the source plane. Similarly, most of the voxels below D c (i.e., point B in Fig. 2) can see the entire detector but only part of the source. To simulate SBDX tomosynthesis reconstruction analytically from the prior CT volumes using the BAA model, we directly apply the blurring to the CT images on a plane-byplane basis and sum the blurred planar data into a single tomosynthesis image as where the two rectangular functions denote the blurring kernel, and a(z,z ′ ) denotes the amplitude function Here, C is a constant that preserves the total signal. We can ignore C because the images will be normalized before the subtraction. d(z,z ′ ) denotes the geometric distortion function The amplitude function is proportional to the inverse of the geometry distortion, because the amplitude function preserves the total amount of the blurring structure when the geometric distortion is applied. Both the amplitude function and the geometry function have two parts because the voxels above and below D c have different x-ray paths. Additionally, each of the rectangular window functions represents the finite sizes of the detector and source planes. The blurring kernel is a function of the plane-to-source distance z and the plane-toplane distance z ′ − z, and the amplitude of the blurring kernel is proportional to the inverse of the plane-to-plane distance. When z = z ′ , the blurring kernel becomes a delta function; thus, there is no blurring in the focal plane. In Fig. 2, the paths of the x-rays that pass through the points A and B are not vertically symmetric. The blurring kernel at the central ray is not equal to the off-center blurring kernel. In other words, the out-of-plane blurring is shiftvariant. Therefore, we apply the geometric distortion function in Eq. (12) to modify the shape of the blurring kernel such that the model fits the true blur over the entire plane. After applying the distortion function to each plane pair, we can use convolution or FFT to compute the BAA-simulated image as the sum of the blurred slices where h(x,z ′ ;z) is the blurring kernel of the plane at z ′ relative to the tomosynthesis plane at z. The blurring kernel can be calculated using the term inside the square brackets of Eq. (10) or using tomographic angles along the central ray.
The tomographic angle φ(x,z) is given by And the approximated blurring kernel along the central line is In the 3D case, the 2D blurring kernel of the SBDX system is separable in the x-z and y-z dimensions Combining all components, the BAA model is illustrated in Fig. 3. The inverse geometry tomosynthesis system results in a series of 2D rectangular blurring kernels. In the BAA F. 3. Schematic of the proposed fast convolution-based BAA model with a rectangular blurring kernel for the SBDX system. Each plane is first geometrically distorted following the center lines determined by the system geometry and plane z location. Geometrically distorted planes are then convolved with an appropriately normalized blurring kernel determined by the system geometry and source-detector pair trajectory. The convolved plane data are summed in the z direction to simulate tomosynthesis reconstruction at a specified focal plane. model, Eq. (12) is an ideal geometric distortion function, which ensures that the blurred voxel shadow falls on the correct locations along the x-ray paths. Each pair of planes in the geometric distortion function requires at least one 2D image interpolation. Thus, in Eq. (12), an image volume with n z planes requires O(n 2 z ) 2D image interpolations in total. However, this large number of interpolations can be reduced by interpolating all slices to the same scale of the plane at D c using the distortion functioñ Then, one can compute the convolution and summation as described in the BAA model. The last step involves reinterpolation of the planes back to their original scales using the inverse of the distortion functiond(z,z ′ ). This computation only requires O(n z ) 2D interpolations, which is a significant reduction in computational complexity and time.

2.D. OPAST
The goal of this study is to use a prior CT volume to suppress out-of-plane structures that are the main artifacts in the SAA reconstruction images of the SBDX system. Because the BAA model describes tomosynthesis images on a planeto-plane basis, it is also capable of computing the undesired out-of-plane artifacts in the SAA reconstructed plane. Assuming that we have accurate 3D registration results between the SAA reconstructed image and the BAA simulated image, the out-of-plane artifacts can be suppressed by subtracting simulated artifacts from the reconstruction where f art (x, y,z;t) is the simulated artifact using the blurand-add model. The subtraction technique is greatly affected by mismatched image intensities. Hence, we normalize the image volumes to zero mean and unit variance; µ recon and σ recon denote the mean and standard deviation of the SAA images in Eq. (8), respectively; µ sim and σ sim denote the mean and standard deviation of the BAA simulated images in Eq. (13), respectively. The BAA-simulated out-of-plane artifacts are given by where t denotes the number of neighboring planes that will be preserved during the subtraction. The out-of-plane structures caused by the neighboring planes within the preserved thickness t∆ z will not be subtracted from the reconstructed images f recon . The slice profile of the OPAST image is sharper if a smaller thickness is selected to be preserved. Moreover, instead of using a sharp step-like function to distinguish between the preserved structures and out-of-plane artifacts, we suggest using a smooth function to avoid rapid changes in the depth profile, for example, where the value of k controls how much out-of-plane artifact will be subtracted. Figure 4 presents the depth profiles of round objects in the SAA reconstruction and OPAST images with different k values. A smaller k value leads to sharper depth profiles of the objects. The OPAST algorithm has more obvious effects on large objects than on small objects. A limit of the OPAST algorithm is that the subtraction images are sensitive to the mismatch between the reconstruction and the simulated artifacts. There are three main sources of the artifacts in the OPAST image: error in the BAA model, mismatches caused by patient motion, and noise in the SAA reconstruction. Using a smaller value of k reduces the magnitude of remaining signals in the OPAST image. However, because the magnitude of the out-of-plane artifacts increases as the slices distance decreases, the magnitudes of the errors caused by the BAA model and patient motion increase. The prior CT volume and the BAA simulated artifact images are assumed to be noiseless and independent to the quantum noise in the SAA images. Thus the magnitude of the quantum noise in the SAA reconstruction is not affected by the subtraction. Although small k values can provide sharp depth profiles, the relative magnitude of the artifacts to the remaining true structures may be too large. Therefore, the value of k should be selected to balance the newly introduced artifacts and the out-of-plane artifacts.
Moreover, the shape and the position of the nodule are likely to be different between the reconstruction and the BAA simulated images due to patient motion. When the simulated artifacts include nodule information, subtraction will affect the visual shape and position of the nodule, especially in the z direction. One approach to eliminate the effect of the subtraction artifacts on the nodule is to remove the target nodule from the prior CT image and then simulate nodule-free out-of-plane artifacts. Using this approach, the out-of-plane artifacts caused by other chest structures can still be suppressed without affecting the visual position of the nodule. The nodule removal approach can be used when the motion artifact is very strong.

2.E. Fast 3D image registration
The real-time 3D-3D image registration is an important step to ensure proper execution of the OPAST algorithm. The BAA model can simulate SBDX images with extended FOV from the prior CT volumes. Then the registration becomes the task of finding the window position in the BAA-simulated image that best matches the real-time reconstructed image assuming that the image registration with respect to patient position and rotations is performed once at the beginning of the procedure under human supervision. The automatic 3D-3D image registration is performed only for shifts in the x, y, and z directions in real-time.
Correlation coefficients (CC) are used as a measure of alignment for the image registration. This measurement of similarity was computed between the reconstruction and the windowed BAA-simulated images. The CCs are computed on a plane-by-plane basis using 2D FFT, and then the average To obtain a clear maximum of the CC in the z direction, we propose using the subtraction method to improve the depth resolution during the image registration. The basic idea of subtraction-based image registration is: with a rough estimation of alignment in the z direction, a small number of out-of-plane artifacts from distant planes can be safely subtracted (using a large value of k). Then, the algorithm provides a more confident estimate of the alignment in the z direction (a sharper peak in the CC measurements) using the OPAST images. In the second and subsequent iterations, the out-of-plane artifacts from closer planes can then be subtracted. Thus, the alignment confidence is increased with increasing iterations. The following code is the pseudocode of the subtraction-based iterative image registration using different k values: • Compute CC between f recon (x,y,z) and f sim (x,y,z).
• Find best alignment (x shift (z),y shift (z)) for each z shift .
• Find best temporary vertical alignment coordinatesž shift .
As shown in Fig. 5(b), reducing the k values creates sharper CC curves over the shifts in the z direction. To clarify, the registration in the x and y directions is usually adequate using the unsubtracted images. Hence, the CC over the x and y directions could only be calculated once at the beginning.
The rotation and deformation of the patient are currently ignored in the real-time registration process. A more complex registration method that is able to capture patient deformation and motion can possibly reduce mismatches between the simulated BAA images and the SBDX reconstructions in the OPAST method. However, adding rotation in the prior CT volume requires generating a new set of BAA images.

SIMULATIONS
We performed numerical simulations to evaluate the BAA model and the OPAST algorithm for the SBDX system. The first set of experiments was designed to evaluate the accuracy of the BAA simulation model through comparison with a simulated reconstruction using projections and backprojections. The shift-variant effect was also evaluated by adding displacement to the CT volumes. The second set aimed to evaluate the effects of the OPAST images in the ideal scenario, i.e., no shift-variant effect, motion, or deformation. The third set was designed to evaluate the OPAST and image registration algorithms with simulated respiratory motions.
We used three 3D thorax CT volumes as test cases in the numerical simulations. All three CT volumes were acquired from in vivo lung cancer patients and included pulmonary nodules with 10-15 mm diameters. The parameters of the CT datasets are listed in Table I. The first test case has a well-circumscribed nodule near the diaphragm. The nodule in T I. The parameters of CT datasets used in the numerical simulations. The third CT dataset is upsampled by 2 in the axial direction to achieve a more uniform resolution. the second test case is directly connected to the pleural surface. The pulmonary nodule in the third test case is connected to vessels near the heart. The CT images were first shifted such that target nodules were at the center of the FOV. The volume size of the SBDX reconstruction images was 128 × 128 × 32 with a spacing of 0.5 × 0.5 × 3 mm in the x, y, and z directions, respectively. The BAA simulated image used a volume size of 128 × 128 × 52 with the same spacing as the reconstructions. The average detector signal was approximately 60 counts/pixel, which is typical for SBDX imaging at 5 frames/s. A Poisson random variable was used to simulate the projection data, with no additive noise (e.g., electronic). The effective x-ray energies were 50 keV for both the SBDX and the CT systems.
We used the CC, mean squared error (MSE), and structural similarity index 28 (SSIM) to assess the accuracy of the BAA analytical model and the image quality of the OPAST images. The CC and MSE were computed for the entire 3D volumes after being normalized to a mean of zero and unit standard deviation. The SSIM combines three similarity components: luminance, contrast, and structure. We scaled the images to a mean of 128 and standard deviation of 32 and used the same parameters as those presented in Ref. 28. We calculated the average SSIM value over all pairs of planes to represent the SSIM of the 3D volume. The CT images were downsampled to 3 mm slice thickness to match the tomosynthesis reconstructions and OPATS images in the comparison.

3.A. BAA model
To demonstrate the accuracy of the BAA model, we compared images using both the SAA reconstructions from simulated SBDX projections and the direct BAA simulations using the same CT data. The SBDX is a shift-variant system in which voxels at different positions have different blurring kernels. Therefore, the accuracy of the BAA model will decrease as the center of the reconstructed volume moves away from the center of the simulation volume. To evaluate the errors caused by the shift-variant effect, we simulated reconstruction images with shifts of the location of the nodule relative to the center of the BAA model from −25 to +25 mm in the x, y, and z directions separately. The CC, MSE, and SSIM were then calculated between the reconstructed images and the corresponding volumes in the BAA simulated images with perfect registration.

3.B. OPAST algorithm
To test the proposed OPAST algorithm, we generated a set of OPAST images with different values of the parameters k in Eq. (20). We first aligned the BAA-simulated images with the SAA reconstructions to ensure that the artifacts were subtracted from their correct positions. In addition, we also simulated the OPAST images with the nodule removed (OPAST-NR) as described in Sec. 2.D. The iterative morphological filtering method was used to segment out and remove the nodules from the CT data. 29 Subtraction of the noduleremoved artifacts only suppresses the out-of-plane artifacts caused by the structures other than the nodules. Therefore, the nodule image, including the depth profile, should not be affected by the subtraction.
The reconstruction and OPAST images were compared with the original CT images. However, the SAA reconstruction is the simplest reconstruction method used in digital tomosynthesis. Out-of-plane structures are blurred out only due to the tomographic angle without any additional filtering or processing. We also investigated the maximum likelihood (ML) method to reconstruct images from the same projection data. This method is based on the realistic assumption that photon counts are Poisson random variables. 30,31 The algorithm suppresses the out-of-plane artifacts using the known projection matrix by performing an iterative matrix inversion.

3.C. OPAST with respiratory motions
It is clear that the OPAST method must work in the presence of respiratory motions. The patient motions generate mismatches between the true object and the prior CT image. These mismatches can cause artifacts in the OPAST images. Furthermore, the image registration algorithm must also be robust enough to find a good alignment for the subtraction. To form a sequence of moving images, we warped the CT volumes according to the motion functions in Zeng et al. 16 for a 4-s breathing cycle. The maximum motion in the entire CT volume is 3-4 cm, and the maximum motion of the nodules in the datasets is about 2 cm. The amounts of added nodule motions can be visualized in Fig. 10.
The respiratory motions were applied to the reference CT images to generate new projections and SAA reconstructions. To test the performance of the OPAST method in the presence of respiratory motion, we used the BAA model to generate simulated images as well as artifact images from the CT images without added motion. The registration and subtraction were performed between the simulated images without motion and the reconstructed images with motion. Again, we used the CC, MSE, and SSIM to evaluate the OPAST results. The OPAST images were compared with the CT images with motion. We performed the evaluation for a sequence of motions with a 0.2-s time step for a 4-s breathing cycle. Figure 6 presents the accuracy assessments of the BAA model using the three patient CT datasets. Because the BAA model is shift-variant, the accuracy of the BAA model decreases as the centers of the SAA reconstruction volume and the BAA simulation volumes move away from each other. The shift in the z direction has less effect on the accuracy than those in the x and y directions because shifting in the z direction has a small effect on the tomographic angles and on the relative geometric distortions among planes. For all three similarity metrics, the displacement between the centers of the two volumes does not greatly affect the accuracy of F. 6. Image similarity metrics between BAA-generated and SAA-generated tomosynthesis volumes calculated as a function of shift in the SAA model away from the central ray of the BAA model. Accuracy assessments of the BAA model and shift-variant effects are shown in the x, y, and z directions. The registration was known perfectly when calculating the similarity measurements. The first row uses CT dataset 1, the second row uses CT dataset 2, and the third row uses CT dataset 3. the BAA model for shifts ranging from −10 to +10 mm in both the x and y directions. The CC values within this region are larger than 0.99, the MSEs are smaller than 0.02, and the SSIMs are larger than 0.998. Therefore, we set ±10 mm as the range in which we can ignore the shift-variant effect for the BAA model. Displacements larger than this range may require either moving the patient or resimulating the BAA image with the lesion correctly positioned at the center of the field of view. This range can also be interpreted as a 20 × 20 mm 2 tolerable region of the BAA model. In later results, we will demonstrate that errors caused by the BAA model are not major sources of error in the OPAST image quality. respectively. In the coronal view, the subtraction images provide clearer shapes of the pulmonary nodules than the alternative reconstructed images. The OPAST images with a smaller suppression parameter k provide improved visualization of the nodule and other lung structure details.

4.B. OPAST algorithm
The quantum noise in the reconstructions becomes more apparent if more out-of-plane artifacts are subtracted from the reconstruction (OPAST-2). In the axial view, the vertical resolution has been significantly improved by the subtraction. Column (g) shows the OPAST-NR image with k = 4 and the nodule removed from the simulated out-of-plane artifacts. The nodules become more apparent in both the coronal and axial views; however, the depth resolution of the nodule is similar to that of the SAA reconstructions. The last column (h) displays simulated projections in which the nodules are difficult to identify. Table II presents the results of the image quality evaluations of the reconstruction and OPAST images for three test cases. In general, the OPAST images have much better similarity measurement results than the SAA and the ML reconstructions for all three test cases. Smaller k values yield better agreement between the OPAST images and the CT data. The OPAST-4 images are slightly better than the OPAST-2 images in most cases because the quantum noise and OPAST artifacts become (a) (   Figure 8 presents the subtraction results corresponding to maximum respiratory motion in coronal and axial views. The OPAST images in columns (c)-(e) can track motion of the nodule and lung structures that are showed in column (a). The OPAST images with patient motion can provide similar nodule visibility as the OPAST images without patient motion in columns (f)-(h). The boundaries of the nodules are more blurred in the OPAST images with motion but still visible.

4.C. OPAST with respiratory motions
There are also more structures in columns (c)-(e), which may be caused by the mismatch and motions. Figure 9 shows the reconstruction, subtraction, and artifact images on the same scale. For smaller k values, the magnitude (standard deviation) of the OPAST images is smaller, and the magnitudes of the subtracted out-of-plane artifacts and the motion artifact are larger. The errors caused by the BAA model (including the shift-variant effects) are at the edges of the anatomic structures. The errors caused by the respiratory motions are much stronger than the BAA model errors and quantum noise. Most of the errors caused by the motion artifacts are at a low frequency that may not greatly affect the nodule visibility. Figure 10 presents three image quality measurements of the OPAST images over a breathing cycle for the three test cases. The image quality decreases when the large motions are added to the reconstructions. In some cases, the OPAST-2 images are even worse than the SAA images because the mismatches caused by the motions generate very strong artifacts. In addition, when the motions are large (t near 2 s), the OPAST-8 yields better quality measurements than the OPAST-4. The OPAST-8 is more robust to the respiratory motions. In some extreme cases, the nodule removal method can be applied. The similarity measurements also fluctuate because the motions lead to different structures and motion artifacts in the OPAST images. Figure 10 column (d) shows the added respiratory motions (dotted lines) of the lung nodule and the measured motions (solid lines) of the entire volume using the proposed image registration method. The registration method can accurately track the motions in the x and y directions. In the z direction, the registration method is limited by the volume spacing and the small tomographic angle of the system.

DISCUSSION
In the current OPAST algorithm, the entire reconstructed region is treated as a single static volume during the registration. In other words, the deformations within the reconstruction volume were ignored when compared with the BAA simulated image. However, the motion-affected CT images were generated using interpolation such that the sizes and shapes of the bones, nodules, vessels, etc., were very different from the reference, which causes more difficulties in the image registration and subtraction. Moreover, the BAA-simulated image and artifacts were also assumed to be constant during the procedure and were not recalculated unless the real reconstruction volume was determined to be significantly shifted from the center of the simulation. Inhale and exhale motions of the lung were considered shifts of the entire reconstruction volume. The motions of the structures outside of the FOV may also affect the results of the registration. The closer voxels are likely to have less relative motions to the volume than the distant voxels. According to Eq. (15), the magnitude of the out-of-plane artifacts is proportional to the inverse of the planar distance. Therefore, the structures from planes that were distant from the VOI have less of an effect than the planes within or next to the VOI. More robust algorithms such as using deformation maps and mutual information could further benefit the OPAST method although computation time would increase. The BAA subtraction model can be generalized to any tomosynthesis system with application of geometry specific derivations. The proposed methods are particularly useful on the SBDX system because the SBDX system contains a very large number of projections (2500) and operates in real-time. The OPAST method uses a simple motion model and registration algorithm, and therefore the image quality achieved may not surpass other correction algorithms with CT priors. However, the algorithm is well-suited to the real-time requirement of the SBDX system since the computational demand is low and the required image processing can keep up with image acquisition. With some modifications, the subtraction method can be potentially used in other applications in which real-time image guidance is required.

CONCLUSION
In this paper, we describe a method for out-of-plane artifacts suppression using the BAA model for the SBDX system. The BAA model was derived from the projection-tobackprojection physical model and used to simulate images and out-of-plane artifacts from a patient-specific prior CT image. The accuracy of the BAA model and the shift-variant effect were evaluated using numerical simulations. The nodule visibility was improved by the OPAST method, especially in the z direction. The image registration and subtraction were also tested with added respiratory motions. Our proposed OPAST algorithm offers the advantages of suppressing artifacts using previously acquired CT data while also providing the near-real-time temporal resolution required for TBNbx image guidance using the SBDX system.