Frontiers in the Convergence of Bioscience and Information Technology
View this Special IssueResearch Article  Open Access
Multimodal Data Integration for ComputerAided Ablation of Atrial Fibrillation
Abstract
Imageguided percutaneous interventions have successfully replaced invasive surgical methods in some cardiologic practice, where the use of 3Dreconstructed cardiac images, generated by magnetic resonance imaging (MRI) and computed tomography (CT), plays an important role. To conduct computeraided catheter ablation of atrial fibrillation accurately, multimodal information integration with electroanatomic mapping (EAM) data and MRI/CT images is considered in this work. Specifically, we propose a variational formulation for surface reconstruction and incorporate the prior shape knowledge, which results in a level set method. The proposed method enables simultaneous reconstruction and registration under nonrigid deformation. Promising experimental results show the potential of the proposed approach.
1. Introduction
Current treatment of cardiac arrhythmias ranges from noninvasive strategies, such as pharmacological therapy, to minimally invasive techniques, such as catheterbased ablation, and to open surgical techniques. While medical therapy can mitigate the occurrence of arrhythmias, these treatments may have significant side effects since most drugs used have some toxicity that is not suitable for longterm therapy. The catheterbased procedure is proven to be an effective method in treating patients with certain cardiac arrhythmias [1]. It is much less invasive and more established. It also demands shorter recovery time than the surgical approach. Thus, catheterbased radio frequency (RF) ablation has become a widely accepted method in the treatment of cardiac arrhythmias, including atrial fibrillation (AF) and ventricular tachycardia (VT). These arrhythmias affect a large number of people and result in significant morbidity and mortality.
AF is the most common sustained cardiac arrhythmia encountered in clinical practice. In the United States alone, there are over 3.5 million patients with this disorder [2]. AF can result in serious complications, including congestive heart failure and thromboembolism. Despite recent advances, drug therapy to control this disease is still unsatisfactory. As an alternative, a nonpharmacological, interventional approach based on creating percutaneous catheterbased lesion inside the heart has been developed. Lesions are delivered in the left atrialpulmonary vein junction with an aim to electrically isolate these veins from the rest of the atrium. This protects the atrium from fast heart beating that is originated in the veins, which initiate and perpetuate AF.
The procedure of interventional AF treatment entails mapping the left atrium and the attached pulmonary veins using an electroanatomic mapping (EAM) system. This mapping information can be used to deliver lesions as well. This electrical approach is suitable for the heart since it is an electromechanical organ, where mechanical contractions are driven by electrical stimulus. However, there is a serious limitation of the EAM system in that it is not able to provide an accurate anatomical information of heart. Typically, a virtual shell is used to represent the atrial wall and the vein. The points on the atrial wall, where the catheter is manually touched, are used to create this shell [3].
The catheterbased ablation process can be greatly improved if a real anatomy is used instead of the virtual shell. To ensure safe catheter maneuverability and enable delivery of effective lesions with minimal collateral damage and complications, it is critical to have both the anatomical information and the electrical information available to the operator. This is particularly important for performing ablation in a complex structure such as the left atrium that is surrounded by important organs, which are vulnerable to damage if lesions are not appropriately directed with close anatomical guidance. Furthermore, even the pulmonary veins themselves are liable to be damaged with grave longterm consequences if the lesions extend deeply into the veins instead of being restricted to the ostia.
The use of a multimodal data integration process can provide an anatomical, physiological, and functional representation simultaneously. In practice, this can be achieved by combining an anatomical surface model acquired by MRI/CT images and the localized electrical information measured by an EAM system. In the registration process, one obvious difficulty stems from the noise and/or outliers that are inevitably associated with the MRI/CT imaging process and the EAM data collection procedure. Unlike other organs in the body, heart undergoes contractile motion, apart from respiratory motion, thus making it unique and very challenging to register and integrate data of different modalities. In addition to physiological variations such as changes in the heart rate, the heart rhythm, and the respiratory effect, various types of heart motion are the source of outliers.
The main contribution of this work is to provide enhanced imaging of the anatomical heart surface from sparse and noisy EAM data in combination with a heartshape model obtained from MRI/CT reconstruction as a prior knowledge. For 3D surface reconstruction, we propose a variational formulation in the level set framework that is an efficient numerical scheme. The level set method is particularly of great use in representing a shape due to its topologyfree and implicit characteristics. By leveraging the 3D heartshape model, we can compensate incomplete EAM data, thereby representing the anatomical heart more accurately. The proposed method has two important advantages. First, it is robust against nonrigid deformation caused by cardiac motion and noise. Second, it can construct the optimal surface without an explicit correspondence between the MRI/CT surface and EAM data due to the implicit surface representation.
The rest of this paper is organized as follows. Previous related work is reviewed in Section 2 followed by the proposed multimodal data integration method for computeraided ablation of atrial fibrillation in Section 3. Both synthetic and real data are tested to demonstrate the efficiency of the proposed method in Section 4. Concluding remarks and future work are given in Section 5.
2. Review of Related Work
Developing a computerguided system for ablative heart surgery involves image registration or integration techniques. They are usually performed under a rigid transformation between preoperative MRI/CT reconstruction and intraoperative EAM data points [4, 5]. Among various registration algorithms, the iterative closest point (ICP) method and its variants have been widely used for this application due to their computational efficiency [6].
The ICP algorithm begins with two meshes and an initial guess for their relative rigidbody transform. It refines the transform iteratively by generating pairs of corresponding points on the meshes and minimizing an error metric repeatedly [7]. However, the standard ICP algorithm does not take noise and outliers into account. Since noise and outliers may affect the ICP performance substantially, several ICP variants have been proposed in [8] to mitigate this problem. One popular approach to identify outliers is to use a threshold, including a certain constant, a fraction of a sorted distance, and some multiple of the standard deviation of a distance [9–11]. Even with these variants, it is still challenging to deal with nonrigid deformation and differentiate inliers from outliers.
Most of previous schemes used an ICPbased method without addressing the abovementioned problem. Instead, they focused on clinical registration. For example, Reddy et al. [5] and Malchano [6] showed the feasibility of combining MRI with CARTOXP in a porcine model of myocardial infarction (MI). They used the modified Iterative Closest Point (mICP) scheme for registration, but did not address the outlier problem. The modification is to adopt hierarchical registration by adding the class information in the algorithm. A clinical registration strategy that combines landmarks and surface registration was proposed in [4]. This study assessed the accuracy for each cardiac chamber using a different clinical registration method. It was observed in [12] that the size of the left atrium affects the accuracy. The patient who has a bigger chamber volume tends to have more ablation errors.
The rigid transformation assumption made by existing schemes is simple yet insufficient in most cases. It often yields unsatisfactory results since a nonrigid deformation is involved between the anatomical heart model reconstructed by MRI/CT images and temporal instances of the heart at the collection of EAM data points. This physiological and anatomical variation that occurs in the formation of the heart surface model and the collection of EAM data points demands a nonrigid transformation (or equivalently diffeomorphism) between the model and the data. Woo et al. proposed a novel image integration technique by incorporating nonrigid deformation using the level set method in [13].
To overcome the limitation of the traditional registration approach based on the rigidtransformation assumption, we formulate this problem as a 3D surface reconstruction problem from EAM data points with a given surface prior. A similar context arises in surface reconstruction from point clouds in a scanned noisy image. Surface reconstruction using an explicit representation has been considered by researchers, for example, [14, 15]. Typically, this approach needs to parameterize a large point set that could be difficult to manipulate. Another approach was proposed in [16, 17] to construct triangulated surfaces using Delaunay triangulations and Voronoi diagrams. It has to determine the right connection among points in the point set, which could be challenging in handling noisy and unorganized point data.
Surface reconstruction based on an implicit shape representation using the level set technique has been studied for almost two decades by applied mathematicians, for example, [18–20]. The nonparametric (or implicit) surface representation has an advantage in dealing with arbitrary topology change and deformation. Hoppe et al. [21] proposed an algorithm in reconstructing a surface using the signed distance function from unorganized points. Zhao et al. [22] proposed another algorithm using the unsigned distance function and the weighted minimal energy to reconstruct the surface. These algorithms are however restricted to situations where the population of points is dense enough to characterize the target surface. They are not applicable to our application where EAM data points are sparse and insufficient.
The problem of insufficient EAM data encountered in the 3D heart surface construction, including the left atrium and its pulmonary veins, can be mitigated by incorporating a heartshape prior provided by MRI/CT imaging. Then, the optimal surface can be obtained by minimizing the energy functional that consists of a datafitting term and a prior knowledge term as detailed in the next section.
3. Multimodal Data Integration for Simultaneous Surface Reconstruction and Registration
In this section, we present a multimodal data integration algorithm for simultaneous surface reconstruction and registration. This algorithm reconstructs the heart surface from measured EAM data points and a heartshape prior obtained by MRI/CT imaging using the level set method.
3.1. Surface Reconstruction and Registration
Under the level set framework [23], a surface in (a curve in ) to reconstruct can be represented by the zerolevel set of a higher dimensional embedding function as given by where is the domain of . This implicit representation can provide the geometric shape of surface with the region defined by the union of and its interior, denoted by . Then, the shape of is given by where is the heaviside function as defined by
For the representation of a surface model obtained from MRI/CT images, we define its shape in a similar way using another level set function as given by . Now, we denote the set of measured EAM data points by where is the number of data points.
The essential assumption in this surface reconstruction application is that surface to reconstruct is an equivalent class to a given prior surface model under small nonrigid deformation and the surface is close to the data points in . Thus, surface can be obtained by minimizing an energy functional that consists of a data fitting term and a prior knowledge term. Our goal is to find the embedding function associated with surface that minimizes the following energy functional: where detail of each term will be described in the following section.
3.2. Derivation of Energy Terms
The energy functional in (4) consists of two terms. The first term measures how well the surface is fit to measured points based on the distance between the surface and these points. The second term measures how plausible the surface is in terms of the prior knowledge of the target surface. Parameter is a weight that adjusts the importance of these two factors.
The data fitting term can be written as where is Dirac measure and this term sums up the Euclidean distance between point and . Recall that is a signed distance function where the value of each point gives the euclidean distance between the point and the interface. To impose the prior knowledge on the target surface, should be close to prior surface model with a smoothsurface assumption. In other words, we penalize any abrupt change of the surface gradient. Here, we use two terms to represent the prior knowledge, that is, where is the smoothness regularization term in the form of and is the shape dissimilarity term of the following form: and where is a rigid transformation resulting from scaling, rotation, and translation. Note that the smoothness regularization term measures the length of the 2D curve and the area of a 3D surface while the shape dissimilarity term measures the symmetric difference between and under a rigid transformation.
Combining (4)–(8), the total energy functional is given by Then, surface reconstruction can be achieved using the energy minimization principle as under constraint , which is a property of a signed distance function. Instead of applying the same weight for both and as shown in (9), different weights can be used for each term.
3.3. Numerical Implementation
For numerical implementation, we use the following approximations for the heaviside function and the Dirac delta measure as defined in [24, 25]:
The energy functional in (10) can be minimized with respect to using the EulerLagrange equation with defining the initial surface. Finally, the gradient descent method is applied to the resultant EulerLagrange equation, which leads to where denotes the exterior normal to the boundary , and denotes the normal derivative of at the boundary.
To discretize the equation in , we adopt a finite differences explicit scheme. The usual notations are as follows: let be the space step and let be the time step. The finite differences are expressed as
We first set as the initial surface and then update using the following discretization:where
Solving the above partial differential equations numerically is challenging since the time step should be constrained to a small value in maintaining numerical stability. In addition, it is computationally expensive to find a high dimensional surface. Thus, it is desired to employ an efficient numerical scheme and we naturally use multigrid method that adopts a hierarchical representation of the data in multiple scales and propagates the solution from the coarse scale to the fine scale to achieve computational efficiency.
3.4. Relationship Between Variational Formulation and Bayesian Inference
This variational approach presented in Sections 3.1 and 3.2 can be interpreted from the interpretation of Bayesian inference under a probabilistic framework. This relationship is presented in this subsection. The target surface can be obtained by maximizing the following posterior probability: where is the likelihood function, is the prior probability of the surface. Maximizing this conditional probability with data point for surface is equivalent to minimizing its negative logarithm where is a constant. Thus, by setting we can convert (17) to (4).
4. Results and Discussion
We begin with simple yet illustrative examples to demonstrate the efficiency and robustness of the proposed algorithm. We use synthetic geometric objects in 2D and 3D which have geometric features that aim to be preserved under reconstruction process. Then, the evaluation of the algorithm is performed based on real patient data.
4.1. Synthetic Data
We first compare the proposed scheme with the ICP scheme in registration accuracy using a synthetic 2D star shape as shown in Figure 1. The original 2D star shape (image size: pixels) is shown in Figure 1(a). It is deformed as shown in Figure 1(e). Furthermore, 33 noisy contour points are generated by adding Gaussian noise (2%, 6%, and 10% standard deviation of contour points, resp.) to original points obtained by sampling the original shape in Figure 1(a). The deformed shape stands for the reconstructed heart shape and data points with different noise levels represent the EAM data points. EAM data points can have errors from sensor, heart movement, and patient breathing. In synthetic experiments, we used Gaussian noise to represent noises introduced in EAM data points.
(a) Original shape
(b) ICP (2% noise)
(c) ICP (6% noise)
(d) ICP (10% noise)
(e) Deformed shape
(f) Proposed (2% noise)
(g) Proposed (6% noise)
(h) Proposed (10% noise)
Graphical illustration of the registration results between deformed shape of different noise level and noisy contour points are presented in Figures 1(b)–1(d) for the ICP scheme. In Figures 1(f)–1(h), shape reconstruction using both shape prior Figure 1(e) and noisy contour points are presented using the proposed scheme. It is clear from Figure 2 that the proposed scheme outperforms ICP. For quantitative error analysis, we measure the mean Euclidean distance between the reconstructed surface and measured data points with varying degree of Gaussian noise. The result is shown in Figure 2. Again, the proposed algorithm outperforms ICP significantly. This is especially true when the noise level is higher. The proposed algorithm produces more stable mean Euclidean distance than ICP as well.
Next, we compare the proposed scheme with ICP using a 3D synthetic jar example. Experimental results are shown in Figure 3, where the original and the deformed shapes are shown in Figures 3(a) and 3(f), respectively. Points extracted from the corrupted surfaces with various Gaussian noise levels (3%, 6%, 9%, and 12%) are used for visual evaluation. Registration between deformed surface and data points with different noise levels using ICP is shown in Figures 3(b)–3(e). Reconstructed surfaces based on data points at different noise levels with the proposed algorithm are presented in Figures 3(g)–3(j). The proposed method generated the reconstructed surface which incorporates the data point as well as deformed prior shape in Figure 3(e). The mean distances are also measured for accuracy comparison as shown in Figure 4. Again, the proposed algorithm is significantly better than the ICP scheme in terms of stability and distance, which is especially obvious at higher noise levels, as shown in Figure 4.
(a) Original shape
(b) ICP (3% noise)
(c) ICP (6% noise)
(d) ICP (9% noise)
(e) ICP (12% noise)
(f) Deformed shape
(g) Proposed (3% noise)
(h) Proposed (6% noise)
(i) Proposed (9% noise)
(j) Proposed (12% noise)
We can adjust the importance of both a data fitting term and a prior knowledge term that includes a regularization and shape similarity term by tuning the parameter . For the choice of the parameter , we use that gives reasonable results. In general, if is too small, then the importance of the data fitting term increases, which results in overfitting to the data point. On the other hand, too big produces the reconstructed shape which is almost same as the prior shape. Thus, parameter is desired to be carefully chosen according to the application.
4.2. Patient Data
The final example is a set of real patient data. 3D preoperative contrastenhanced MR angiography (MRA) was performed to delineate endocardial boundaries of the left atrium and pulmonary veins. The voxel size was mm and 45 slices were used in the experiment. We obtained MRA and 250 EAM data points from the same patient. The EAM data consists of the CARTO points imported from the CARTOXP, including measurement points as well as ablation points.
After delineating and removing unwanted regions such as the left ventricle (LV) and other small veins, we reconstruct the 3D model as shown in Figure 5 using ITKSNAP [26] and Matlab software. Afterwards, a twostep registration process is applied for the ICP scheme. First, we perform the landmark registration using three junctions between LA and pulmonary veins: LALIPV, LALSPV, and LARSPV. These points are used for the initial pose of subsequent registration. Second, surface registration using the ICP scheme is performed to refine accuracy furthermore. The resulting image is shown in Figure 6(a).
(a) MRA of LA
(b) 3D Reconstruction
(a) ICP
(b) Proposed method
To validate the proposed algorithm, the optimal surface is reconstructed using 250 EAM data points by incorporating a heart shape prior from preoperative MRA. By minimizing the energy functional, the final result is shown in Figure 6(b), where diamonds (in blue) represent EAM data points, and circles (in red) represent ablation points. Blurred points are located inside. A quantitative evaluation result can be obtained in terms of the mean Euclidean distances of EAM and ablation points from the surface of the left atrium. They are reported in Table 1, which shows that the proposed approach outperforms the ICP method significantly.

5. Conclusion and Future Work
A novel multimodal data integration technique using the level set method for catheter ablation of AF was presented in this paper. This technique enables reconstruction and registration simultaneously using data fitting, regularization, and shape prior energy terms. It provides better performance than the existing ICP method in accuracy. In the proposed framework, the heartshape model from MRA reconstruction is used as a prior shape knowledge. Thus, we can use the shape information to compensate for insufficient EAM data. Clinically, this technique can improve efficacy and safety of AF ablation by integrating EAM data and 3D imaging data.
Dynamic cardiac shape analysis will make the current integration method more precise and meaningful. We plan to incorporate a richer set of spatiotemporal shape models using dynamic shape information in the future. Besides, we may consider a localized regularization method around the point data to obtain more precise reconstruction.
References
 D. P. Zipes and H. J. J. Wellens, “What have we learned about cardiac arrhythmias?,” Circulation, vol. 102, no. 4, pp. 52–57, 2000. View at: Google Scholar
 I. B. Ray and E. Heist, “Treating atrial fibrillation: what is the consensus now?,” Postgraduate Medicine, vol. 118, no. 4, pp. 47–58, 2005. View at: Google Scholar
 C. Pappone, S. Rosanio, and G. Oreto et al., “Circumferential radiofrequency ablation of pulmonary vein ostia: a new anatomic approach for curing atrial fibrillation,” Circulation, vol. 102, no. 21, pp. 2619–2628, 2000. View at: Google Scholar
 J. Dong, H. Calkins, and S. B. Solomon et al., “Integrated electroanatomic mapping with threedimensional computed tomographic images for realtime guided ablations,” Circulation, vol. 113, no. 2, pp. 186–194, 2006. View at: Publisher Site  Google Scholar
 V. Y. Reddy, Z. J. Malchano, and G. Holmvang et al., “Integration of cardiac magnetic resonance imaging with threedimensional electroanatomic mapping to guide left ventricular catheter manipulation: feasibility in a porcine model of healed myocardial infarction,” Journal of the American College of Cardiology, vol. 44, no. 11, pp. 2202–2213, 2004. View at: Publisher Site  Google Scholar
 Z. Malchano, “Image guidance in cardiac electrophysiology,” Massachusette Institute of Technology, Boston, Mass, USA, June 2006. View at: Google Scholar
 P. J. Besl and H. D. McKay, “A method for registration of 3D shapes,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 14, no. 2, pp. 239–256, 1992. View at: Publisher Site  Google Scholar
 S. Rusinkiewicz and M. Levoy, “Efficient variants of the ICP algorithm,” in Proceedings of the 3rd International Conference on 3D Digital Imaging and Modeling (3DIM '01), pp. 145–152, Quebec, Canada, MayJune 2001. View at: Google Scholar
 D. Chetverikov, D. Svirko, D. Stepanov, and P. Krsek, “The trimmed iterative closest point algorithm,” in Proceedings of the 16th International Conference on Pattern Recognition (ICPR '02), vol. 3, pp. 545–548, Quebec, Canada, August 2002. View at: Publisher Site  Google Scholar
 T. Masuda, K. Sakaue, and N. Yokoya, “Registration and integration of multiple range images for 3d model construction,” in Proceedings of the 13th International Conference on Pattern Recognition (ICPR '96), pp. 879–883, Vienna, Austria, August 1996. View at: Google Scholar
 K. Pulli, “Multiview registration for large data sets,” in Proceedings of the 2nd International Conference on 3D Digital Imaging and Modeling (3DIM '99), pp. 160–168, Ottawa, Canada, October 1999. View at: Google Scholar
 E. K. Heist, J. Chevalier, and G. Holmvang et al., “Factors affecting error in integration of electroanatomic mapping with CT and MR imaging during catheter ablation of atrial fibrillation,” Journal of Interventional Cardiac Electrophysiology, vol. 17, no. 1, pp. 21–27, 2006. View at: Publisher Site  Google Scholar
 J. Woo, B.W Hong, S. Kumar, I. B. Ray, and C.C J. Kuo, “Joint reconstruction and registration using level sets: application to the computeraided ablation of atrial fibrillation,” in Proceedings of International Conference Frontiers in the Convergence of Bioscience and Information Technologies, Jeju, Korea, October 2007. View at: Google Scholar
 G. C.H. Chuang and C.C. J. Kuo, “Wavelet descriptor of planar curves: theory and applications,” IEEE Transactions on Image Processing, vol. 5, no. 1, pp. 56–70, 1996. View at: Publisher Site  Google Scholar
 D. Rogers, An Introduction to NURBS: With Historic Perspective, Morgan Kaufmann Publishers, San Francisco, Calif, USA, 2001. View at: Google Scholar
 N. Amenta and M. Bern, “Surface reconstruction by Voronoi filtering,” Discrete and Computational Geometry, vol. 22, no. 4, pp. 481–504, 1999. View at: Publisher Site  Google Scholar  MathSciNet
 H. Edelsbrunner, “Shape reconstruction with delaunay complex,” in Proceedings of the 3rd Latin American Symposium on Theoretical Informatics (LATIN '98), pp. 119–132, Springer, Campinas, Brazil, April 1998. View at: Google Scholar
 A. Dervieux and F. Thomasset, A Finite Element Method for the Simulation of a RayleighTaylor Instability, vol. 771 of Lecture Notes in Mathematics, Springer, Berlin, Germany, 1979. View at: Google Scholar
 A. Dervieux and F. Thomasset, “Multifluid incompressible flows by a finite element method,” in Seventh International Conference on Numerical Methods in Fluid Dynamics, W. Reynolds and R. MacCormack, Eds., vol. 141 of Lecture Notes in Physics, pp. 158–163, Springer, Berlin, Germany, 1980. View at: Google Scholar
 S. Osher and J. A. Sethian, “Fronts propagating with curvature dependent speed: algorithms based on the HamiltonJacobi formulation,” Journal of Computational Physics, vol. 79, no. 1, pp. 12–49, 1988. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, and W. Stuetzle, “Surface reconstruction from unorganized points,” Computer Graphics, vol. 26, no. 2, pp. 71–78, 1992. View at: Publisher Site  Google Scholar
 H.K. Zhao, S. Osher, and R. Fedkiw, “Fast surface reconstruction using the level set method,” in Proceedings of the 1st IEEE Workshop on Variational and Level Set Methods, The 8th IEEE International Conference on Computer Vision, pp. 194–202, Vancouver, BC, Canada, July 2001. View at: Publisher Site  Google Scholar
 J. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Cambridge University Press, Cambridge, UK, 1998. View at: Google Scholar
 T. F. Chan and L. A. Vese, “Active contours without edges,” IEEE Transactions on Image Processing, vol. 10, no. 2, pp. 266–277, 2001. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 H.K. Zhao, T. Chan, B. Merriman, and S. Osher, “A variational level set approach to multiphase motion,” Journal of Computational Physics, vol. 127, no. 1, pp. 179–195, 1996. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 P. A. Yushkevich, J. Piven, H. Cody, S. Ho, J. C. Gee, and G. Gerig, “Userguided level set segmentation of anatomical structures with ITKSNAP,” Insight Jounral, vol. 1, 2005, special issue on ISC/NAMIC/MICCAI Workshop on OpenSource Software. View at: Google Scholar
Copyright
Copyright © 2008 Jonghye Woo et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.