Int J CARS (2009) 4:189–202 DOI 10.1007/s11548-009-0283-4
ORIGINAL ARTICLE
Vascular tree reconstruction with discrete tomography: intensity based camera correction for 3D reconstruction C. Bodensteiner · C. Darolti · A. Schweikard
Received: 10 January 2008 / Accepted: 29 December 2008 / Published online: 4 February 2009 © CARS 2009
Abstract Purpose This paper is concerned with the reconstruction of vascular trees from few projections using discrete tomography. However, its computational cost is high and it lacks robustness when the data are inconsistent. We improve robustness by incorporating an intensity-based cameracorrection method. The proposed approach is also capable of handling small motion artifacts by modeling them as repositionings of a virtual X-ray camera. We also present a parallel implementation which substantially reduces reconstruction time. Methods We propose a data-driven reduction of positional inconsistencies by minimizing the reconstruction residual to increase the robustness. Inspired by motion compensation algorithms in SPECT imaging, we combine an intensity-based 2D/3D-registration method with iterative reconstruction methods. Our objective is the robust vascular-tree reconstruction from positionally inconsistent data. The speed of the reconstruction is substantially increased by a volume-splitting scheme that allows parallel processing. Results Vascular trees in the liver can be accurately reconstructed from few positionally inconsistent projections using digitally reconstructed radiographs. We have tested the proposed method on synthetic projection data and on objects imaged with a new robotized C-arm. We measured a decrease in the average reconstruction residual of about 13% for real data compared to projection data without preprocessing. Over 4,600 reconstruction experiments were conducted to evaluate the speed-up obtained when employing the volume-splitting scheme. Reconstruction time decreased linearly with C. Bodensteiner (B) · C. Darolti · A. Schweikard Institute of Robotics and Cognitive Systems, University of Luebeck, Luebeck, Germany e-mail:
[email protected]
increased number of processor-cores, both for real and synthetic data. Conclusions The proposed method reduces inconsistencies caused by positioning errors and small motion artifacts. No prior segmentation or detection of correspondences between projections is necessary, because all algorithms are intensitybased. As a result, the proposed method allows for robust, high-quality reconstructions, while reducing radiation dose substantially. Keywords 3D Reconstruction · Digital subtraction angiography · Algebraic reconstruction · 2D/3D-registration · Camera calibration
Introduction Three dimensional reconstruction of vascular trees from few projections has been of great interest to researchers over the past decades [4,12,26]. It remains an area of active research, as reasons include lack of robustness and prohibitive runtime of conventional methods. Reconstruction methods can be roughly divided into geometry- and model-based methods and intensity-based methods. Geometry and model-based methods [4,26] rely on matching features in multiple 2D X-ray projection images by imposing geometrical constraints. This involves the non trivial tasks of medical image segmentation [12] and finding of correspondences. An advantage of these techniques is that the reconstruction and camera calibration tasks [5,15] can be combined. However, ambiguous correspondences at vessel crossings and bifurcations in the projections pose major challenges. Intensity-based methods do not need any finding of correspondences. We focus on the intensity-based discrete
123
190
tomography (DT) technique [11], which intrinsically models the discrete nature of vascular tree reconstruction. The reconstruction object is modeled as a function of three spatial variables. The function value denotes the absorption coefficient at a 3D point in the imaged object. In case of DT, this function must be discrete, i.e., its range must consist of a few discrete values. In medical imaging, a discrete-valued function can be obtained with digital subtraction angiography (DSA). For DSA, pairs of projections are acquired before (mask image) and after injection of contrast agent. After subtracting the mask image from the contrast-agent image, the resulting projection data ideally show only the X-ray absorption caused by the contrast agent in the vessels. This a-priori information enables high-quality 3D-reconstructions from very few projections. Additionally, the projection data can be taken from a limited angle range, which is a highly desirable property for intra-operative imaging. Intensity-based reconstruction methods are also affected by a number of problems. Mobile devices (C-arms) are prone to imaging errors due to inherent positioning inaccuracies caused, for example, by mechanical flex. These slight mismatches between the recorded and the actual position of the imaging device have a negative impact on 3D reconstruction results. Compared to reconstructions from fixed devices, the reconstruction errors from mobile devices are much larger due to positional inconsistencies. Other sources of inconsistencies are, for example, image noise and patient motion. In this paper, we propose a data-driven reduction of positional inconsistencies, in order to increase the robustness of 3D-reconstruction from images acquired using mobile devices. The proposed approach is based on an intensitybased camera-correction method. Inspired by motion compensation algorithms in SPECT imaging [13,21], we combine a rigid 2D/3D-registration algorithm with iterative reconstruction algorithms to create a method for simultaneous intensity-based vascular tree reconstruction and correction of positional inconsistencies. This approach is capable of correcting small motion artifacts by modeling them as the repositioning of a virtual X-ray camera. We plan to use this reconstruction technique in combination with a robotized C-arm for intra-operative imaging and navigation. Specifically, we aim at soft tissue navigation applications which use vascular trees for the registration of pre-operative planning data. The structure of this paper is as follows. First, we present the algorithms involved in inconsistency reduction and 3D-reconstruction, and show how they are employed in the designed method. Second, we provide the details of a novel volume-splitting scheme for an efficient parallelization of the reconstruction. Finally, we validate the proposed methods with over 4,600 numerical experiments on simulated as well as real projection data.
123
Int J CARS (2009) 4:189–202
Methods Algebraic 3D-reconstruction methods have the property of calculating an average solution of the reconstruction problem (1). We use this property to correct for positional inconsistencies of the X-ray camera. The camera’s position is given by its extrinsic parameters. We iteratively approximate these parameters using an intensity-based rigid 2D-3D registration of the original projections with the former algebraically reconstructed volume. A new algebraic reconstruction is then started with the corrected projections. After convergence in this correction phase, the corrected projections are used for DT reconstruction of the final volume. Figure 1 shows the schematic overview of this method. Subsequently, we describe the steps of the proposed method in detail. DT techniques can be roughly divided into three categories [11]: heuristic approaches, methods based on the discretization of classical algorithms (e.g. binary steering [3]), and algorithms based on optimization. For an excellent survey of DT and its applications in medical imaging we refer to [11]. The DT algorithm employed in this paper is based on optimization. The reconstruction problem is modeled as a linear programming (LP) problem which can easily incorporate prior information, as detailed in the Appendix. We chose to base the implementation of the DT-algorithm on LP, because LP-algorithms are known to be efficient in solving large scale optimization problems. Formally, 3D-reconstruction with algebraic methods can be formulated as the following linear least-squares problem: find x ∈ Rn such that R(x) = Ax − b2 = min!,
(1)
in which A denotes the projection-operator producing cone beam projections, x the reconstruction volume and b the measured projections. Let f i (η; b; ·), F(η; b; ·) : Rn → Rn be the applications of: x, ai − bi ai , i = 1, . . . , m, ai 2 F(η; b; x) = ( f 1 ◦ · · · ◦ f m ) (η; b; x). f i (η; b; x) = x − η
(2) (3)
Here ai denotes the ith row of the system matrix A, bi the ith element of the measured projection data b, and η is a relaxation parameter. Using the algebraic reconstruction technique (ART) algorithm [2] we can now calculate a limit point for the series k = 0, 1, 2, . . .: x k+1 = F(η; b; x k ),
(4)
∈ e.g. = 0. Censor et al. [2] with starting point proved that, for sufficiently small values of η, the algorithm converges to a least-squares solution (1) with respect to the projection equations. Consider now that the X-ray camera transformations can be modeled as 3D rigid object motion Φi ∈ R6 with six x0
Rn ,
x0
Int J CARS (2009) 4:189–202
191
Fig. 1 Schematic overview of the iterative correction process: obtain an initial volume with an algebraic 3D-reconstruction. Iterate correction and algebraic reconstruction. In each iteration, correct the camera posi-
tion using 2D/3D registration, then algebraically reconstruct the volume from the corrected projections. Finally, use the corrected projections for a DT reconstruction
degrees of freedom. In situations of positional errors, we want to minimize the following residual reconstruction error:
are digitally reconstructed radiographs (DRR) generated from the previously reconstructed volume through virtual X-ray cameras at the new positions. This algorithm is then iteratively repeated until convergence. The main building blocks of the 2D/3D-registration algorithm are a distance measure and an optimization scheme [19]. The distance measure is intensity-based and robust since it does not require prior segmentation. Specifically, it is a gradient-based correlation distance measure [19]. The cost function to be optimized is calculated as the average of the normalized cross-correlations of the partial derivatives of the projection and DRR image intensities
J
x Φ
=
M
Ai x( Φi ) − Bi 2 = min!.
(5)
i=1
Here Ai denotes the subset of projection equations that produce one cone beam projection image. Bi denotes the corresponding subset of X-ray measurements and M the number of acquired projection images. The minimization procedure is iterative, alternating camera parameter estimation and algebraic reconstruction. In each iteration, we first estimate the parameters Φi of Eq. (5) using 2D/3D registration of the current reconstruction result with the original projections. We then solve Eq. (1) with the above ART reconstruction algorithm. In the 2D/3D registration step, 3D-transformations of the positions of the X-ray cameras are calculated, such that a distance measure between the acquired projections and a set of forward projections is minimized. The forward projections
C x =
x,y (Pr x (x,
y) − Pr x )(Drr x (x, y) − Drr x ) , 2 2 x,y (Pr x (x, y) − Pr x ) x,y (Drr x (x, y) − Drr x )
(6) C = 21 (C x + C y ), C y analog to C x . Pr y (x, y) = and Pr x (x, y) =
∂ I Pr (x,y) ∂x
∂ I Pr (x,y) ∂y
denote the partial derivatives of
123
192
image intensities I Pr (x, y), and Pr x and Pr y denote their mean values. The optimization is carried out in a hierarchical multi-level approach: a pyramid of images is constructed by down-sampling the image at multiples of two and the optimization is run for each level of the pyramid. This allows for a robust and fast optimization, since evaluating the cost function at lower levels is computationally very cheap. We plotted the values of the cost functions for transformation parameters varying around the optimum found by the 2D/3D registration algorithm. The plots show a relatively smooth, monotonically increasing function with a unique maximum and a good convergence range. They also show that out of plane parameters (e.g. movements in the viewing direction of the X-ray camera) cannot be determined as accurately as the other parameters [25]. As a result, we penalize strong out-of-plane parameter corrections by downscaling the component of the correction vector in the direction of the central ray. In addition, we fix the first acquired image during correction in order to prevent a movement of the reconstruction volume due to the erroneous determination of the out of plane parameters. A variant of the Hooke–Jeeves pattern-search [18] was used for 2D/3D parameter optimization. To avoid double function evaluations the function values were cached in a hash table as described in [7]. Additionally, the camera-correction algorithm has a nice geometric interpretation. Assuming the reconstructed volume has n voxels and a projection image has m pixels, the projection equations can be interpreted geometrically as hyperplanes in an n-dimensional image space. The iterative reconstruction algorithm ART shifts the current solution in the direction of hyperplane normals towards a least squares solution with respect to all hyperplanes defined by all projection images. Subsequently, the 2D/3D registration algorithm searches for a rigid transformation of a subset of m hyperplanes defined by a projection image and its X-ray camera position. The sought transformation minimizes the distance from the current least-squares solution to this subset alone.
Int J CARS (2009) 4:189–202
A new least-squares solution with respect to all transformed hyperplanes is then calculated by the algebraic reconstruction algorithm.
Volume-splitting scheme for parallelized implementation The combined reconstruction and registration method is implemented in C/C++ and parallelized with OpenMP. Since mobile C-arms usually acquire cone beam projection data on a circular trajectory (e.g. the C-arm rotates around the isocenter in a circular orbit), we propose a simple method to enhance parallelization by splitting the reconstruction task into two independently solvable sub-problems. A sub-problem is obtained by partitioning the reconstruction volume using a plane. This splitting plane is defined by the iso-center and the rotation axis (normal vector of the splitting plane) of the C-arm rotation during object imaging. An example of a splitting plane is shown in Fig. 2. The resulting two reconstruction regions can then be solved in parallel. The reconstruction problem will not be partitioned into parts of equal size, since the form of the peel volume determines the unknowns and, therefore, the computational complexity of the reconstruction. The partitioning is dependent on the reconstruction object and not known in advance, but a significant reduction of reconstruction time is achievable even with disproportionate splits. We provide numerical results of the scaling behavior of this parallel method in the “Results” section. All software libraries are implemented or ported to 64bit code to allow memory usages above 4GB, since LP reconstruction needs huge amounts of memory. This property of LP-based reconstruction algorithms is also a limiting factor for volume resolutions above 2563 voxels. The forward projection calculation of the rigid registration preprocessing is based on a shear warp factorization algorithm [25]. This step takes under 35s/15s for a 2563 /1283 volume (3 down sampling levels) with X-rays of size 568 × 568 pixels. For the DT reconstruction we used the barrier algorithm of the multi-threaded LP library CPLEX 10.0.
Fig. 2 Visualization of a splitting plane (denoted in green). Left splitting of the peel volume. Right splitting plane in the resulting reconstruction result
123
Int J CARS (2009) 4:189–202
Experiments This section describes the setup of the experiments conducted with synthetic and real X-ray images. We have used liver data sets from three different patients. As these patients were not imaged with our C-arm, projections were simulated. For this reason we refer to these experiments as synthetic projection data experiments. We include two other clinical data sets—a thorax CT and a distal femur bone—in order to show that the positional error correction algorithm is not restricted to DSA data. The method was also tested on data sets acquired with our prototypical C-arm. The test objects are three complex phantoms and a chess-figure data set, as described in Sect. “Robotic C-arm projection data”. The synthetic and C-arm projection data were employed in over 4,600 reconstruction experiments with various testing goals: reconstruction quality, reconstruction speed and robustness in case of positionally perturbed projection data. We thus provide a short outline of the conducted experiments. First we turn our attention to simulated projection data and analyze reconstruction from few (3–5) projections with respect to error and runtime. These experiments give insight into the properties and the drawbacks of the DT reconstruction technique. We then perturb the extrinsic parameters in order to generate synthetic inconsistent projections. The robustness of the proposed method with respect to inconsistencies is analyzed by reconstructing from the perturbed synthetic projections. This analysis is important since mobile C-arm devices are prone to positional errors causing inconsistent projections. The accuracy of reconstructions for the synthetic data experiments was estimated using the following standard error measures that are defined by voxel differences between a reconstructed volume v and a ground truth volume v [11]: 2 vi − vi 2 vi − vi , VE = SE = vi + vi vi + vi vi − v i MV = . (7) 2 vi The volume error V E measures only the volume difference between the reconstructed volume v and a ground truth volume v. This could be problematic in a scenario in which an object is reconstructed as too thin in some places and too thick in others. The remaining distance measures are normalized versions of the number of misclassified voxels. The misplaced voxels (M V ) measure [20] favors volumes with a low total number of non-zero voxels. The shape error (S E) [17] measure is more sensitive to wrongly reconstructed voxels due to the factor 2 in the nominator. In our last experiments, we test the position correction properties of the proposed method for four different phantoms imaged using
193
our robotized C-arm prototype. For these data sets, we present visual reconstruction results obtained using only a few of the acquired radiographs. The accuracy of the reconstruction cannot be measured with the indicators (7) because the ground truth volume is not available, since we cannot obtain a perfect reconstruction from the C-arm device. The following sections describe the details of the experiments using synthetic and real projection data. Synthetic projection data The simulated radiographs (5682 pixel) were generated from the CT-data of livers from three different patients. The vascular trees were manually segmented by MeVis1 for liver resection planning. Projections were generated from the segmented vascular trees and the original tree was reconstructed from a few projections selected from this data. We refer to the vascular trees in Fig. 3 as A for the left, B for the middle, and C for the right dataset. In general, vessel trees with many vessel crossings/branchings on the projection images lead to higher dimensional—and therefore more difficult—optimization problems. Therefore, vascular tree B is more difficult to reconstruct than vascular tree A. Vascular tree C has a higher spatial resolution (2563 voxels with size 0.72 mm) and poses a larger optimization problem. To create the reconstruction from few projection scenarios, we varied the number of projections used for the reconstruction and the angle between them. For this purpose, we simulated cone beam projection data (40 projections, 180◦ angle range, 4.5◦ step width) for each vascular tree using an alpha-clipping approach [24]. Figure 3 in Sect. “Synthetic projection data” shows the vascular trees and examples of the generated forward projections. Since the outcome of the discrete tomography reconstructions is a volume with values x ∈ [0..1], we must convert it to a binary volume x ∈ {0, 1}. We thus apply a simple threshold (0.5) to the volumes in order to set all voxels that are greater than the threshold to one, and zero otherwise. This subsequently allows the voxelwise comparison of volumes. The parameters for the influence of the regularization terms and the error tolerance were set to α = 0.3, β = 0.2, τ0 = 0.5, τ1 = 0.5 in all experiments. The intrinsic camera parameters were chosen analog to the real mobile X-ray device: focal length 960 mm, pixel scale 0.40 mm/pixel. Robotic C-arm projection data Real projection data were acquired with an experimental robotized C-arm (Ziehm–Imaging). Given the direct and 1
The segmentation was carried out by MeVis (http://www.mevis.de), a company which provides expert segmentations of liver vascular trees for liver resection procedures.
123
194
Int J CARS (2009) 4:189–202
Fig. 3 Vascular trees for simulation reconstruction tests. Left a 1283 voxels with size 1.8mm, middle b 1283 voxels with size 1.8 mm, left c 2563 voxels with size 0.72 mm. Bottom corresponding DRRs, 5682 pixels, pixel scale 0.40 mm
inverse kinematics, the C-arm can be accurately positioned to take images at arbitrary positions [1]. Details about the kinematic solution and the calibration of the robotic C-arm device are given in [14], which are beyond the scope of this paper. We used a frame-grabber device for the image acquisition, which allows for a software-based magnification and rotation of the projection images. Geometrical distortions of the X-ray camera were determined using a grid phantom attached to the image amplifier. The calibration object is shown in Fig. 4. The distortion model is based on bi-variate polynomials of degree five as in [7]. The polynomial coefficients were calculated a priori. Figure 4 also shows images that have been undistorted using the computed coefficients. We used four different test objects for the experiments. Three objects were phantoms constructed by inserting PVC-tubes in Agarose gel, as shown in Fig. 5. The PVC tubes served as a simplified model of a blood vessel system and can be easily filled with contrast agent. The fourth object was a chess piece inserted in a plastic container. For this highly non-convex object, we reconstruct the shape of the contrast agent around it. This type of object poses a very difficult reconstruction scenario; it is therefore a challenge for the proposed method. The known object shape allows for a good visual judgment of the reconstruction quality. We imaged these setups from different positions before and after filling of contrast agent. For every experiment we used 50 ml Imeron 300 (active pharmaceutical ingredient: Iomeprol).
123
Results The results section is structured in simulated and real data, in analogy to the “Experiments” section. Reconstructions with simulated projection data The original volume represents the ground truth for simulated projection experiments. The quality of the 3D-reconstructions can thus be easily evaluated from simple voxel wise comparisons with the standard error measures described in Sect. “Experiments”. The number of projections used for the reconstruction and the angle between them were varied as described in Sect. “Synthetic projection data”. We reconstructed from all projection subsets of size 2 for vascular trees A and B, and plotted the error of every combination in Fig. 6. The figure shows that the reconstruction error drops very fast if the angle between the projections is greater than 15 degrees. The percent of average misplaced voxels was below 3.1%, which represents a very small deviation from the ground truth. A maximal misplaced voxel error of 0.17% was observed for vascular 40 tree A, and 3.1% for vascular tree B, in a number of 2 = 780 reconstruction experiments. The results of these experiments are summarized in Table 4. The reconstruction error increases with the complexity of the shape, as can be seen for tree B. The angle between the projections also determines the size of the peel volume and, therefore, the problem size, i.e., the number of unknowns. This correlation is exemplary plotted
Int J CARS (2009) 4:189–202
195
Fig. 4 Geometrical distortion correction with bivariate polynomials of degree 5. Left original projection (bottom) and distortion corrected projection (top). Middle, right de-warped projections by using the attached calibration grid
Fig. 5 Experimental setup for the robotized C-arm. a A phantom made of Agarose gel with inserted PVC-tubes. b Projection image without contrast agent, c projection image with contrast agent, d corresponding DSA image
Fig. 6 Plot of the misplaced voxels (%) versus the angle between two projections. Every point in the plot corresponds to a reconstruction experiment with two projection images. Left vascular tree A.
Right vascular tree B. The angle between projections is calculated by taking the difference between the projection number times 4.5◦ (|Pr oj Nr.1 − Pr oj Nr.2|∗ 4.5◦ )
in Fig. 7 for vascular tree A, in which problem size (number of unknowns + number of constraints) is plotted versus the angle of two projections. Vascular tree C could be accurately reconstructed with 3–5 projections (see Table 5) as well.
(BIF), Soft Bounds (SB), Regularized Best Inner Fit (RBIF), Regularized Soft Bounds (SB). These algorithms are described in the Appendix. Vascular tree complexity plays a major role in reconstruction quality and computation time, as can be seen by analyzing Tables 1, 2, 3, 4 and 5. The term error voxels in the tables refers to all corresponding voxels between a ground truth volume and a second volume. It denotes the sum of all voxels which are set in the reconstruction volume, but not in
Speed-up of parallel DT reconstruction algorithms The reconstruction time depends on the employed DT-reconstruction algorithm, and increases in the order Best Inner Fit
123
196
Int J CARS (2009) 4:189–202
Fig. 7 Every point in the plot corresponds to a reconstruction problem with two projection images. Left plot of the number of unknowns versus angle range of two projections. Right number of non zero rays versus the angle range of two projections (vascular tree A)
Table 1 Reconstruction results with synthetic projection data and subsets of size 3, vascular tree A and 4 Threads
Table 2 Reconstruction results with synthetic projection data and subsets of size 3, vascular tree B and four Threads
Table 3 Reconstruction results with synthetic projection data and subsets of size 3, vascular tree B and splitted version with 2-2-Threads
Table 4 Reconstruction results with synthetic projection data and subsets of size 2, optimization function SB from Eq. (17), and splitted version with 4-4 Threads)
123
Optimization function
Max/min/med error voxels
MV (%)
SE (%)
VE (%)
Reconstruction time max / min / med
#
BIF
567 / 1 / 35
0.6
1.2
0.5
200 / 15 / 51
179
R-BIF
198 / 1 / 6
0.1
0.2
0.1
900 / 37 / 161
179
SB
4/1/1
0.02
0.04
0.03
262 / 29 / 72
178
R-SB
4/1/1
0.02
0.04
0.03
389 / 58 / 130
178
Optimization function
Max/min/med error voxels
MV (%)
SE (%)
VE (%)
Reconstruction time max / min / med
#
BIF
657 / 154 / 248
1.8
3.5
2.5
471 / 132 / 194
175
R-BIF
375 / 154 / 198
1.4
2.8
2.5
1636 / 391 / 609
175
SB
306 / 181 / 231
1.6
3.3
1.9
1223 / 216 / 347
175
R-SB
286 / 177 / 220
1.6
3.1
2.0
1858 / 417 / 684
175
Optimization function
Max/min/med error voxels
MV (%)
SE (%)
VE (%)
Reconstruction time max / min / med
#
BIF
647 / 154 / 257
1.8
3.7
2.3
428 / 120 / 187
174
R-BIF
370 / 154 / 207
1.5
2.9
2.4
1157 / 289 / 508
174
SB
327 / 181 / 239
1.7
3.4
1.8
957 / 163 / 293
174
R-SB
286 / 177 / 222
1.6
3.2
1.9
1573 / 344 / 538
174
MV (%)
SE (%)
VE (%)
Reconstruction time max / min / med
#
Vascular Tree
Max/min/med error voxels
A
237 / 0 / 10
0.2
0.3
0.1
1428 / 30 / 175
780
B
1424 / 77 / 441
3.1
6.2
1.1
3862/ 161/ 619
780
C
79631/33279/54090
104.1
36.6
73.2
2432 / 95 / 172
317
Int J CARS (2009) 4:189–202 Table 5 Reconstruction results with synthetic projection data for vascular tree C, optimization function SB from Eq. (17) and splitted version with 4-4 Threads
197
Subset size
Max/min/med error voxels
3er subsets
47313/ 0 / 3817
2.6
5.4
4er subsets
17756/ 0 / 795
0.5
1.1
5er subsets
268 / 0 / 18
0.0
0.0
0.0
the ground truth volume and vice versa. The reconstruction time is measured in seconds. The last column (#) denotes the number of experiments. As can be seen in Table 1, we measure an average reconstruction time of 51s with the BIF algorithm, and 72s with the SB algorithm for the simple vascular tree A. Regularized variants needed on average 161s with the R-BIF algorithm and 130s with the R-SB algorithm. The more complex tree B leads to far more difficult and larger problems. The errors are larger than for tree A, and there is a significant increase of reconstruction time, as can be seen in Tables 2 and 3. As a result, the average reconstruction time increased to 194s with the BIF algorithm and 347s with the SB algorithm, respectively, 609s with the R-BIF and 684s with the R-SB regularized variants. This is already too long for intraoperative imaging. However, a very positive result is that all algorithms scaled very well by the use of several processor cores, as shown in Table 1. This observation is also nicely depicted in Fig. 8. We may conclude that a nearly linear improvement in execution time was measured with the multi-threaded versions. Another conclusion is that regularization improved the results qualitatively. However, reconstruction time increases substantially, since the individual reconstruction problems become considerably larger than in the non-regularized case. Every 6-connected (resp. 26-connected) voxel pair in the peel volume introduces a regularization term. A significant reduction of reconstruction time can be as well achieved with the splitting approach for the regularized versions. This is due to the fact that regularization constraints can be well split. This result can be seen in Table 3. For vascular tree C (resolution 2563 , voxel size 0.72mm) we need three projections to achieve very good reconstruction results, since this dataset has a higher spatial resolution. This can be seen in Tables 4 and 5. Table 5 shows that, in this case, five projections are sufficient to achieve perfect reconstructions. Reconstructions from disturbed projection data For the camera position-correction experiments, the camera Cam to T pertCam . positions were artificially perturbed from TWorld World The translation parameters were perturbed by maximally ±7 mm and the rotation parameters were perturbed by
MV (%)
SE (%)
VE (%)
Reconstruction time max / min / med
#
1.5
9865 / 216 / 568
438
0.2
2363 / 512 / 968
76
1688 / 540 / 794
56
maximally ±1◦ . First, we registered 150 positional distorted DRRs to the original volume and calculated the Euler angle representation of the deviation matrix T m to validate the 2D/3D registration algorithm World Cam T m Cam regCam = TregCam TWorld .
(8)
The average registration accuracy showed deviations of 0.137 (x), 0.134 (y) and 0.27 (z) degree and an average translational deviation of 0.25 (x), 0.80 (y), 0.35 (z) mm, i.e. the errors were very small. We subsequently performed the proposed reconstruction algorithm on the perturbed data. For the DT reconstructions we used the error measures previously introduced. We conducted 12 experiments for vascular tree A and nine experiments for vascular tree B. By perturbing the position of the projection data, the volume position and orientation are also slightly modified. This prohibits a direct comparison of reconstructions with the ground truth. We therefore performed a rigid registration of the original volume and the reconstruction volume to compensate for this effect. However, the 3D/3D-registration introduces non-zero voxels at vascular tree boundaries (interpolation and partial volume effects) leading to a distortion of the voxel-wise evaluation. We measured an average misplaced voxel error of 18.8% for vascular tree A and 17.6% for vascular tree B. However, in both cases, the visual results were very good, as can be seen in Figs. 9 and 10. The correction algorithm was also tested on two clinical CT datasets. Figure 11 shows a clear reduction of inconsistencies. Note that the correction method is not restricted to vascular trees and works very well for standard CT datasets as well. We mention that the reconstruction residual decreased by 29.32% for the distal femur dataset and by 13.42% for the thorax CT dataset, relative to the uncorrected reconstructions. Reconstruction from acquired projection data The correction method works very well with real projection data due to the robust gradient correlation distance measure. In general, a significant improvement of the reconstruction quality is visible after 3–7 iterations of the correction algorithm. We used five projection images taken from angles ranging over 100◦ . To appreciate the quality of the proposed
123
198
Int J CARS (2009) 4:189–202
Fig. 8 Reconstruction time speedup through parallelization with synthetic projection data, 50 instances—3er subsets from vascular tree B. Left SB (Eq. 17). Right R-SB (Eq. 16)
Fig. 9 Reconstruction results for a liver vascular tree using five positionally disturbed projection images. From left to right ART reconstruction from perturbed data (a), DT reconstruction from corrected data (b), original dataset (c)
Fig. 10 Reconstruction results for a liver vascular tree using five positionally disturbed images. From left to right ART reconstruction from perturbed data (a), DT reconstruction from corrected data (b), original dataset (c)
reconstruction, observe the improved edges and higher details in the visualization in Fig. 12. The Agarose phantoms were reconstructed from ten equally spaced projections taken from angles ranging over 110deg. Again, the visualization of the reconstruction in Fig. 13 clearly shows improved edges and better details. For the three Agarose datasets we achieved a reduction of the error residuals, over the uncorrected reconstructions, by 13.13, 12.72 and 12.51% after seven correction iterations.
123
Small deviations in the camera positions lead to significant reconstruction artifacts which are clearly reduced by the positional correction algorithm.
Discussion and future work Vessel tree reconstruction from X-ray projections is a difficult problem. Geometry-based approaches which rely on
Int J CARS (2009) 4:189–202
Fig. 11 Reconstruction results for clinical CT datasets using 100 positionally disturbed projection images: Slice of an ART reconstruction from 100 perturbed projection images of a distal femur bone; same slice of the reconstruction result before a and after b the correction
199
algorithm. Reconstruction from 100 perturbed projection images of a thorax CT (only the inner circular region was reconstructed); same slice of the reconstruction result before c and after d the correction algorithm
Fig. 12 Camera correction results for the knight chess piece. The first three images show the results of the correction algorithm after 0, 3, 5 iterations. The right image shows the DT reconstruction result from the camera corrected data using the RSB approach
Fig. 13 Reconstruction results for the Agarose phantom (II) using ten projection images. From left to right. Reconstruction a without error correction from the experimental C-arm robot. Magnification of a struc-
ture that shows strong reconstruction artifacts b. c shows the same structure after the correction algorithm (7 iterations, 12.5% reconstruction residual decrease)
matching features need to first segment the vessels in the projection images. To this date, no solution for accurate automatic vessel segmentation exists [6]. Segmentation errors
significantly hamper the accuracy of the reconstruction achievable with feature-based approaches. The automatic segmentations can be manually improved to substantially
123
200
increase the quality of the information for reconstruction. In addition to being time consuming and intra operatively not feasible, manual improvement still cannot guarantee good reconstructions, since the most important feature-correspondence problem [26] still needs to be solved in order to automatically calculate correct calibration- and correction parameters. We implemented a livewire semi-automatic segmentation method [16], which allows for accurate segmentation of vessels and vessel centerlines in projection images. A vessel crossing was manually identified in one projection. Using epipolar constraints and multiview geometry, as common in geometry-based approaches [22,26], over 50 possible matches have been identified in a second projection. The correspondence problem is highly ambiguous and therefore very difficult to solve. Even when using additional projection images, a satisfactory reconstruction could not be obtained for the complex vascular tree B which has many multiple vessel crossings. In contrast, we have obtained a discrete tomography reconstruction with the proposed method, which is almost indistinguishable from the ground truth. The mathematical foundations of 3D reconstruction by discrete tomography have been thoroughly researched [3,8–11,24,23], but there are still very few practical implementations for clinical use, e.g., [11]. This is due to the fact that discrete-tomography methods are computationally expensive, or that they are not robust enough with respect to positional and motion errors and artifacts. Very large computation times and memory requirements have thus made discrete tomography unsuitable for intra-operative use in the past. The proposed method addresses both the robustness and the efficiency issues, showing that it is possible to correct positional inconsistencies and motion artifacts in an efficient manner. To our best knowledge, no optimized execution times for reconstruction with discrete tomography have been published. Robert et al. [20] reported computation times ranging from 1 to 6 hours of CPU time on a Sun SparcStation 2 using simulated annealing approach. Weber et al. [23] report 4 min and 6 s for the reconstruction of a single 256 × 256 2D-image using a similar LP-based approach. We have shown that it is possible to reduce the computation time significantly, while still obtaining very good 3D-reconstructions using parallel algorithms. We consider this to be a significant contribution of this paper, since in the near future, multi-core architectures will enable further shortening of reconstruction times. Furthermore, the proposed method is not limited to reconstructing vascular trees. It can deal with different objects with complex geometries for which correspondence finding would be difficult. Current work includes development of a GPU accelerated reconstruction-registration framework and the calibration-reconstruction pre-processing. In vivo experiments with
123
Int J CARS (2009) 4:189–202
porcine livers are scheduled to test the algorithms on clinical projection data. The unequal mixing of blood and contrast agent in the vascular system at different places and at different times could also pose an important problem for reconstruction. To study this problem, we plan to build more realistic vascular tree phantoms that simulate blood flow and contrast agent mixing.
Conclusion In this paper, we described a new method for an intensitybased camera calibration algorithm which does not rely on segmentation or the identification of correspondences. In all our experiments, the inconsistencies were clearly reduced resulting in significantly improved reconstructions, both quantitatively and visually. Complex liver vascular trees were accurately reconstructed from two to five positionally disturbed projections (DRRs). We have shown that the method is not limited to vascular tree, but it can reconstruct different objects with complex geometries, while correcting positional inconsistencies and small motion artifacts. A new volume-splitting approach was designed for the purpose of algorithm parallelization. Experiments have shown that the gain in reconstruction time is almost linear in the number of employed processor cores. For this purpose, over 4,600 experiments with real and synthetic data sets were conducted. Reconstruction experiments with real acquired DSA projection data using a new experimental mobile C-arm were also successfully accomplished.
Appendix The reconstruction problem is modeled, analog to algebraic reconstruction, as a system of projection equations. Each X-ray measurement leads to one equation in the system; the projection value bi for X-ray i is a linear combination of the n voxel values x j weighted by coefficients ai, j . Considering the notations A := (ai, j ), x = (x j ), b = (bi ), i = 1, . . . , m, and j = 1, . . . , n, the system of equations can be written in matrix form: Ax = b.
(9)
Fishburn et al. [8] proposed the first DT approach based on LP (10) min 0, x, Ax = b, 0 ≤ xi ≤ 1, ∀i.
x∈ n
(10)
The binary constraints xi ∈ {0, 1} have been relaxed to interval constraints in order to account for measurement errors (Eq. 10) and to enable polynomial time algorithms (interior point methods). Gritzmann [9,10] suggested the
Int J CARS (2009) 4:189–202
201
Best Inner Fit (BIF-Eq. 11) method, that replaces Ax = b with the inequality Ax ≤ b to reconstruct the maximal volume satisfying the projection constraints in order to deal with small inconsistencies. In the following equations, the inner product of the vectors a, b is denoted by a, b and e is defined as the vector e := (1, 1, . . . , 1)T . minn − e, x , Ax ≤ b, 0 ≤ xi ≤ 1, ∀i.
(11)
x∈
Weber et al. [24] extended these approaches with regular ization terms j,k |x j − xk | where the sum of distances of adjacent voxel pairs are minimized to favor spatially coherent reconstructions (RBIF-Eq. 12) α x j − xk , Ax ≤ b, 0 ≤ xi ≤ 1, ∀i. minn −e, x+ x∈ 2 j,k
(12) Further enhancements include a binarization functional x, e − x to achieve binary solutions without a rounding scheme. Also the projection equations are supplemented with error variables γi ∈ , λi , γ = (γi ), i = 1, . . . , m to account for noise and X-ray measurement errors (Eq. 13). τ γ if γi ≥ 0, τ0 > 0, τ > 0 (13) λi := 0 i τ1 γi else
respective without regularization as Soft Bounds (SB) method. min n β
x∈[0,1]
m
λi .
(17)
i=1
To reduce the dimension of the optimization it is common practice to constrain the number of unknowns. This is achieved by calculating a peel volume [24] that is the reconstructed volume after removing all voxels which lie on a ray with a zero (resp. below a threshold) measurement, i.e. remove all x j with nj=1 ai, j x j = bi = 0, i = 1, . . . , m. However, if the projection data contains positioning or motion inconsistencies, this reduction is problematic and introduces severe reconstruction artifacts. Therefore it is very important to remove positional inconsistencies prior to this reduction step.
µ 2
By assigning different weights τ0 and τ1 to the error variables, one can also choose between a best inner or a best outer fit solution. In this paper, the following LP will be solved (Eq. 14): min n
x∈[0,1] ,z i,k
m α µ z i,k + x, e − x + β λi 2 2 j,k
⎛ a11 x ⎜ .. ˜ ˜ A = b, A := ⎝ . γ am1
(14)
i=1
⎞ · · · a1n 1 . .. .. ⎟ . .. . ⎠ · · · am,n 1
(15)
subject to 0 ≤ xi ≤ 1, γi ∈ λi ≥ τ0 γi , λi ≥ τ1 γi z i, j ≥ xi − x j , z i, j ≥ x j − xi In this formulation, the regularization term j,k |x j − xk | is replaced by additional slack variables z j,k . Because the created LP problems have a very large scale in the clinical situation, we decided against an automatic determination of a binarisation threshold. We have used the following optimization functions that we denote as regularized Soft Bounds (R-SB) min n
x∈[0,1] ,z i,k
m α z i, j + β λi , 2 i, j
i=1
(16)
References 1. Binder N, Bodensteiner C, Matthaeus L et al (2006) Image guided positioning for an interactive c-arm fluoroscope. In: Computer Assisted Radiology and Surgery (CARS) 2. Censor Y, Gordon D, Gordon R (1983) Strong underrelaxation in kaczmarz’s method for inconsistent systems. Numer Math 41:83– 92 3. Censor Y, Matej S (1999) Binary steering of non-binary iterative algorithms. In: Discrete tomography: foundations, algorithms and applications. Birkhauser, Boston, pp 285–296 4. Chen S, Carroll JP (2000) Three dimensional reconstruction of coronary arterial tree to optimize angiographic visualization. IEEE Trans Med Imag 19:318–336 5. Chen S, Metz CE (1997) Improved determination of biplane imaging geometry from two projection images and its application to three-dimensional reconstruction of coronary arterial trees. Med Phys 24:633–654 6. Condurache AP (2008) Cardiovascular biomedical image analysis. Ph.D. Thesis, University of Luebeck 7. Doetter M (2006) Flouroskopiebasierte navigation zur intraoperativen unterstuetzung orthopaedischer eingriffe. Ph.D. Thesis, Technische Universitaet Muenchen 8. Fishburn P, Schwander P, Shepp L, Vanderbei R (1997) The discrete radon transform and its approximate inversion via linear programming. Discr Appl Math 75:39–61 9. Gardner R, Gritzmann P (1997) Discrete tomography: determination of finite sets by X-rays. Trans Am Math Soc 349:2271–2295 10. Gritzmann P, de Vries S, Wiegelmann M (2000) Approximating binary images from discrete X-rays. SIAM J Optim II(2):522–546 11. Herman G, Kuba A (2003) Medical applications of discrete tomography. In: Proceedings of the IEEE, vol 91 12. Kirbas C, Quek F (2004) A review of vessel extraction techniques and algorithms. ACM Comput Surv 36:81–121 13. Lee KJ et al (1998) Use of forward projection to correct patient motion during spect imaging. Phys Medi Biol 43:171–187 14. Matthaeus L, Binder N, Bodensteiner C et al (2007) Closed form inverse kinematic solution for fluoroscopic c-arms. Advanced Robotics 21 15. Metz CE, Fencil LE (1989) Determination of three-dimensional structure in biplane radiography without prior knowledge of the relationship between the two views. Med Phys 16:45–51
123
202 16. Mortensen EN, Barrett WA (1998) Interactive segmentation with intelligent scissors. Graph Models Image Process 60:349–384 17. Onnasch DGW, Prause GPM (1999) Heart chamber reconstruction from biplane angiography. In: Discrete tomography: foundations, algorithms and applications. Birkhauser, Boston, pp 385–403 18. Parker D (1999) An empirical investigation of the global behavior of several pattern search algorithms. In: Technical Report, Department of Computer Science—University of North Carolina 19. Penney G, Weese J, Little JA et al (1998) A comparison of similarity measures for use in 2d-3d medical image registration. In: MICCAI, p 1153 20. Robert N, Peyrin F, Yaffe MJ (1994) Binary vascular reconstruction from a limited number of cone beam projections. Med Phys 21(12):1839–1851 21. Schumacher H, Fischer B (2007) A new approach for motion correction in spect imaging. Bildverarbeitung fuer die Medizin 22. Vaillant AR, Blondel C, Mal G, Ayache N (2004) Reconstruction of coronary arteries from one rotational X-ray. In: Projection Sequence, Research report, no 5214. INRIA (2004)
123
Int J CARS (2009) 4:189–202 23. Weber S, Schnoerr C, Schuele T, Hornegger J (2006) Binary tomography by iterating linear programs. In: Geometric properties for incomplete data. Springer, Netherlands, pp 183–197 24. Weber S, Schuele T, Hornegger J et al (2004) Binary tomography by iterating linear programs from noisy projections. In: Proceedings of International Workshop on Combinatorial Image Analysis (IWCIA). Springer, Auckland 25. Weese J, Goecke R, Penney G, Desmedt P, Buzug T, Schumann H (1999) Fast voxel-based 2d/3d registration algorithm using a rendering method based on the shear-warp factorization. In: Proceedings of the SPIE, pp 802–810 26. Windyga P, Garreau M, Shah M et al (1998) Three dimensional reconstruction of the coronary arteries using a priori knowledge. Med Biol Eng Comput 36:158–164