# Multimodal Data Integration for Computer-Aided  Ablation of Atrial Fibrillation

Authors: Jonghye Woo; Byung-Woo Hong; Sunil Kumar; Indranill Basu Ray; C.-C. Jay Kuo
Journal: Journal of Biomedicine and Biotechnology
 (2008)

---

## 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.


---

## Body

## 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 Section2 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 S(x) in ℝ3 (a curve in ℝ2) to
reconstruct can be represented by the zero-level set of a higher dimensional
embedding function ϕ(x):Ω→ℝ as given by(1)S={x∈Ω∣ϕ(x)=0},interior(S)={x∈Ω∣ϕ(x)>0},exterior(S)={x∈Ω∣ϕ(x)<0},
where Ω is the domain
of ϕ(x). This implicit representation can provide the
geometric shape of surface with the region defined by the union of S and its interior,
denoted by S¯. Then, the shape of S is given by S¯(x)=H(ϕ(x)), where H(x) is the
heaviside function as defined by
(2)H(x)={1,x≥0,0,x<0.For the representation of a surface modelM(x) obtained from
MRI/CT images, we define its shape M¯(x) in a similar
way using another level set function psi(x) as given by H(ψ(x)). Now, we denote the set of measured EAM data points
by(3)D={p1,p2,…,pn}⊂ℝ3,
where n is the number
of data points.The essential assumption in this surface
reconstruction application is that surfaceS to reconstruct
is an equivalent class to a given prior surface model M under small
nonrigid deformation and the surface is close to the data points in D. Thus, surface S 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 S that minimizes
the following energy functional:(4)E(ϕ)=Epoint(ϕ,D)+αEprior(ϕ),
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 α≥0 is a weight
that adjusts the importance of these two factors.The data fitting term can be written as(5)Epoint(ϕ|D)=∑i=1n∫Ω|ϕ(x)⋅δ(x−pi)|2dx,
where δ(x)=(d/dx)H(x) is Dirac
measure and this term sums up the Euclidean distance between point pi 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, S should be close
to prior surface model M 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,(6)Eprior(ϕ)=Ereg(ϕ)+Eshape(ϕ∣ψ),
where Ereg(ϕ) is the
smoothness regularization term in the form of(7)Ereg(ϕ)=∫Ω|∇H(ϕ(x))|dx,
and Eshape(ϕ∣ψ) is the shape
dissimilarity term of the following form:(8)Eshape(ϕ∣ψ)=∫Ω|H(ϕ(x))−H(ψ(T(x)))|2dx,
and where T(x) 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 H(ϕ) and H(ψ) under a rigid
transformation.Combining (4)–(8), the total energy functional E(ϕ) is given by(9)E(ϕ)=Epoint(ϕ∣D)+α(Ereg(ϕ)+Eshape(ϕ∣ψ)).
Then, surface reconstruction can be achieved using the energy minimization principle as(10)ϕ*=arg⁡min⁡ϕ⁢E(ϕ)
under constraint |∇ϕ|=1, which is a property of a signed distance function.
Instead of applying the same weight α for both Ereg and Eshape 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]:(11)δ(z)={0if|z|>ε,12ε[1+cos⁡(πzε)]if|z|≤ε,H(z)={1ifz>ε,0ifz<−ε,12[1+zε+1πsin⁡(πzε)]if|z|≤ε,The energy functional in (10) can be minimized with
respect to ϕ(x) using the Euler-Lagrange
equation with ϕ(t=0,x)=ϕ0(x) defining the
initial surface. Finally, the gradient descent method is applied to the
resultant Euler-Lagrange equation, which leads to
(12)∂ϕ∂t=−2∑i=1nϕ(x)δ(x−pi)+α[δ(ϕ)div(∇ϕ|∇ϕ|)−2∫Ω(H(ϕ)−H(ψ))δ(ϕ)dx],ϕ(t=0,x)=ϕ0(x)inΩ,δ(ϕ)|∇ϕ|∂ϕ∂n→=0on∂Ω,
where n→ denotes the
exterior normal to the boundary ∂Ω, and ∂Ω/∂n→ 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 h be the space
step and let Δt be the time step. The finite differences are expressed as(13)Δxϕ=ϕi+1,j,k−ϕi−1,j,k2Δx,Δyϕ=ϕi,j+1,k−ϕi,j−1,k2Δy,Δzϕ=ϕi,j,k+1−ϕi,j,k−12Δz,Δxxϕ=ϕi+1,j,k+ϕi−1,j,k−ϕi,j,kΔx2,Δyyϕ=ϕi,j+1,k+ϕi,j−1,k−ϕi,j,kΔy2,Δzzϕ=ϕi,j,k+1+ϕi,j,k−1−ϕi,j,kΔz2,Δxyϕ=ϕi+1,j+1,k−ϕi+1,j−1,k−ϕi−1,j+1,k+ϕi−1,j−1,k4ΔxΔy,Δxzϕ=ϕi+1,j,k+1−ϕi+1,j,k−1−ϕi−1,j,k+1+ϕi−1,j,k−14ΔxΔz,Δyzϕ=ϕi,j+1,k+1−ϕi,j+1,k−1−ϕi,j−1,k+1+ϕi,j−1,k−14ΔyΔz.We first setϕn as the initial
surface and then update ϕn+1 using
the following discretization:(14)ϕi,j,kn+1−ϕi,j,knΔt=−2∑m=1nϕi,j,knδ(x−pm)+αh2δ(ϕi,j,kn)κ−2∫Ω(H(ϕi,j,kn)−H(ψ))δ(ϕi,j,kn)dx,where(15)κ=(Δxϕ2Δyyϕ−2ΔxϕΔyϕΔxyϕ+Δyϕ2Δxxϕ+Δxϕ2Δzzϕ−2ΔxϕΔzϕΔxzϕ+Δzϕ2Δxxϕ+Δyϕ2Δzzϕ−2ΔyϕΔzϕΔyzϕ).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 Sections3.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 S can
be obtained by maximizing the following posterior
probability:(16)P(S∣D)=P(D∣S)P(S)P(D),
where P(D∣S) is the
likelihood function, P(S) is the prior
probability of the surface. Maximizing this conditional probability with data
point D for surface S is equivalent
to minimizing its negative logarithm(17)−log⁡(P(S∣D))=−log⁡(P(D∣S))−log⁡(P(S))+c,
where c is a constant.
Thus, by setting(18)Epoint(ϕ,D)=−log⁡(P(D∣S)),αEprior(ϕ)=−log⁡(P(S)),
we can convert (17) to (4).


## 3.1. Surface Reconstruction and Registration

Under the level set framework [23], a surface S(x) in ℝ3 (a curve in ℝ2) to
reconstruct can be represented by the zero-level set of a higher dimensional
embedding function ϕ(x):Ω→ℝ as given by(1)S={x∈Ω∣ϕ(x)=0},interior(S)={x∈Ω∣ϕ(x)>0},exterior(S)={x∈Ω∣ϕ(x)<0},
where Ω is the domain
of ϕ(x). This implicit representation can provide the
geometric shape of surface with the region defined by the union of S and its interior,
denoted by S¯. Then, the shape of S is given by S¯(x)=H(ϕ(x)), where H(x) is the
heaviside function as defined by
(2)H(x)={1,x≥0,0,x<0.For the representation of a surface modelM(x) obtained from
MRI/CT images, we define its shape M¯(x) in a similar
way using another level set function psi(x) as given by H(ψ(x)). Now, we denote the set of measured EAM data points
by(3)D={p1,p2,…,pn}⊂ℝ3,
where n is the number
of data points.The essential assumption in this surface
reconstruction application is that surfaceS to reconstruct
is an equivalent class to a given prior surface model M under small
nonrigid deformation and the surface is close to the data points in D. Thus, surface S 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 S that minimizes
the following energy functional:(4)E(ϕ)=Epoint(ϕ,D)+αEprior(ϕ),
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 α≥0 is a weight
that adjusts the importance of these two factors.The data fitting term can be written as(5)Epoint(ϕ|D)=∑i=1n∫Ω|ϕ(x)⋅δ(x−pi)|2dx,
where δ(x)=(d/dx)H(x) is Dirac
measure and this term sums up the Euclidean distance between point pi 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, S should be close
to prior surface model M 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,(6)Eprior(ϕ)=Ereg(ϕ)+Eshape(ϕ∣ψ),
where Ereg(ϕ) is the
smoothness regularization term in the form of(7)Ereg(ϕ)=∫Ω|∇H(ϕ(x))|dx,
and Eshape(ϕ∣ψ) is the shape
dissimilarity term of the following form:(8)Eshape(ϕ∣ψ)=∫Ω|H(ϕ(x))−H(ψ(T(x)))|2dx,
and where T(x) 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 H(ϕ) and H(ψ) under a rigid
transformation.Combining (4)–(8), the total energy functional E(ϕ) is given by(9)E(ϕ)=Epoint(ϕ∣D)+α(Ereg(ϕ)+Eshape(ϕ∣ψ)).
Then, surface reconstruction can be achieved using the energy minimization principle as(10)ϕ*=arg⁡min⁡ϕ⁢E(ϕ)
under constraint |∇ϕ|=1, which is a property of a signed distance function.
Instead of applying the same weight α for both Ereg and Eshape 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]:(11)δ(z)={0if|z|>ε,12ε[1+cos⁡(πzε)]if|z|≤ε,H(z)={1ifz>ε,0ifz<−ε,12[1+zε+1πsin⁡(πzε)]if|z|≤ε,The energy functional in (10) can be minimized with
respect to ϕ(x) using the Euler-Lagrange
equation with ϕ(t=0,x)=ϕ0(x) defining the
initial surface. Finally, the gradient descent method is applied to the
resultant Euler-Lagrange equation, which leads to
(12)∂ϕ∂t=−2∑i=1nϕ(x)δ(x−pi)+α[δ(ϕ)div(∇ϕ|∇ϕ|)−2∫Ω(H(ϕ)−H(ψ))δ(ϕ)dx],ϕ(t=0,x)=ϕ0(x)inΩ,δ(ϕ)|∇ϕ|∂ϕ∂n→=0on∂Ω,
where n→ denotes the
exterior normal to the boundary ∂Ω, and ∂Ω/∂n→ 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 h be the space
step and let Δt be the time step. The finite differences are expressed as(13)Δxϕ=ϕi+1,j,k−ϕi−1,j,k2Δx,Δyϕ=ϕi,j+1,k−ϕi,j−1,k2Δy,Δzϕ=ϕi,j,k+1−ϕi,j,k−12Δz,Δxxϕ=ϕi+1,j,k+ϕi−1,j,k−ϕi,j,kΔx2,Δyyϕ=ϕi,j+1,k+ϕi,j−1,k−ϕi,j,kΔy2,Δzzϕ=ϕi,j,k+1+ϕi,j,k−1−ϕi,j,kΔz2,Δxyϕ=ϕi+1,j+1,k−ϕi+1,j−1,k−ϕi−1,j+1,k+ϕi−1,j−1,k4ΔxΔy,Δxzϕ=ϕi+1,j,k+1−ϕi+1,j,k−1−ϕi−1,j,k+1+ϕi−1,j,k−14ΔxΔz,Δyzϕ=ϕi,j+1,k+1−ϕi,j+1,k−1−ϕi,j−1,k+1+ϕi,j−1,k−14ΔyΔz.We first setϕn as the initial
surface and then update ϕn+1 using
the following discretization:(14)ϕi,j,kn+1−ϕi,j,knΔt=−2∑m=1nϕi,j,knδ(x−pm)+αh2δ(ϕi,j,kn)κ−2∫Ω(H(ϕi,j,kn)−H(ψ))δ(ϕi,j,kn)dx,where(15)κ=(Δxϕ2Δyyϕ−2ΔxϕΔyϕΔxyϕ+Δyϕ2Δxxϕ+Δxϕ2Δzzϕ−2ΔxϕΔzϕΔxzϕ+Δzϕ2Δxxϕ+Δyϕ2Δzzϕ−2ΔyϕΔzϕΔyzϕ).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 Sections3.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 S can
be obtained by maximizing the following posterior
probability:(16)P(S∣D)=P(D∣S)P(S)P(D),
where P(D∣S) is the
likelihood function, P(S) is the prior
probability of the surface. Maximizing this conditional probability with data
point D for surface S is equivalent
to minimizing its negative logarithm(17)−log⁡(P(S∣D))=−log⁡(P(D∣S))−log⁡(P(S))+c,
where c is a constant.
Thus, by setting(18)Epoint(ϕ,D)=−log⁡(P(D∣S)),αEprior(ϕ)=−log⁡(P(S)),
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
Figure1. The original 2D star shape (image size: 200×200 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.The shape reconstruction results for a synthetic 2D star
shape, (a) the original shape, (b)–(d) ICP results using different Gaussian
noise levels, (e) the deformed shape, and (f)–(h) results of the proposed
method using different Gaussian noise levels.
(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 Figures1(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.Figure 2
Comparison of mean distances between the reconstructed
surface and measured data points for the 2D star example at different noise
levels using ICP and the proposed algorithm.Next, we compare the proposed scheme with ICP using a
3D synthetic jar example. Experimental results are shown in Figure3, 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.The shape reconstruction results for a synthetic 3D
image: (a) the original shape, (b)–(e) ICP results using different Gaussian
noise levels, (f) the deformed shape, and (g)–(j) results of the proposed
method using different Gaussian noise levels.
(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)Figure 4
Comparison of mean distances between the reconstructed
surface and measured data points for the 3D jar example at different noise
levels using ICP and the proposed algorithm.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 α=1/3 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
was0.78125×0.78125×1.5 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 Figure5 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).The 3D patient data: (a) MRA of LV and (b) 3D
reconstruction result of the given MRA.
(a)
MRA of LA(b)
3D ReconstructionThe surface registration results for the patient data.
(a)
ICP(b)
Proposed methodTo 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 Figure6(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.Table 1
Performance comparison of ICP and the proposed method.
ICPProposedEAM point mean distance4.5087 mm2.4113 mmAblation point mean distance3.2046 mm2.0921 mm


## 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
Figure1. The original 2D star shape (image size: 200×200 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.The shape reconstruction results for a synthetic 2D star
shape, (a) the original shape, (b)–(d) ICP results using different Gaussian
noise levels, (e) the deformed shape, and (f)–(h) results of the proposed
method using different Gaussian noise levels.
(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 Figures1(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.Figure 2
Comparison of mean distances between the reconstructed
surface and measured data points for the 2D star example at different noise
levels using ICP and the proposed algorithm.Next, we compare the proposed scheme with ICP using a
3D synthetic jar example. Experimental results are shown in Figure3, 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.The shape reconstruction results for a synthetic 3D
image: (a) the original shape, (b)–(e) ICP results using different Gaussian
noise levels, (f) the deformed shape, and (g)–(j) results of the proposed
method using different Gaussian noise levels.
(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)Figure 4
Comparison of mean distances between the reconstructed
surface and measured data points for the 3D jar example at different noise
levels using ICP and the proposed algorithm.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 α=1/3 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
was0.78125×0.78125×1.5 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 Figure5 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).The 3D patient data: (a) MRA of LV and (b) 3D
reconstruction result of the given MRA.
(a)
MRA of LA(b)
3D ReconstructionThe surface registration results for the patient data.
(a)
ICP(b)
Proposed methodTo 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 Figure6(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.Table 1
Performance comparison of ICP and the proposed method.
ICPProposedEAM point mean distance4.5087 mm2.4113 mmAblation point mean distance3.2046 mm2.0921 mm


## 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.



---

*Source: 681303-2008-02-14.xml*