Abstract
Image-guided percutaneous interventions have successfully replaced invasive surgical methods in some cardiologic practice, where the use of 3D-reconstructed cardiac images, generated by magnetic resonance imaging (MRI) and computed tomography (CT), plays an important role. To conduct computer-aided 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 catheter-based 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 long-term therapy. The catheter-based 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, catheter-based 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 catheter-based lesion inside the heart has been developed. Lesions are delivered in the left atrial-pulmonary 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 catheter-based 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 long-term 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 heart-shape 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 topology-free and implicit characteristics. By leveraging the 3D heart-shape 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 computer-aided 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 computer-guided 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 rigid-body 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 ICP-based method without addressing the above-mentioned problem. Instead, they focused on clinical registration. For example, Reddy et al. [5] and Malchano [6] showed the feasibility of combining MRI with CARTO-XP 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 rigid-transformation 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 heart-shape prior provided by MRI/CT imaging. Then, the optimal surface can be obtained by minimizing the energy functional that consists of a data-fitting 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 heart-shape 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 zero-level 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 smooth-surface 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 Euler-Lagrange equation with defining the initial surface. Finally, the gradient descent method is applied to the resultant Euler-Lagrange 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 contrast-enhanced 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 CARTO-XP, 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 ITK-SNAP [26] and Matlab software. Afterwards, a two-step registration process is applied for the ICP scheme. First, we perform the landmark registration using three junctions between LA and pulmonary veins: LA-LIPV, LA-LSPV, and LA-RSPV. 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 heart-shape 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.