François Chung
email: francois.chung@inria.fr

Hervé Delingette

Multimodal Prior Appearance Models based on Regional Clustering of Intensity Profiles

Model-based image segmentation requires prior information about the appearance of a structure in the image. Instead of relying on Principal Component Analysis such as in Statistical Appearance Models, we propose a method based on a regional clustering of intensity profiles that does not rely on an accurate pointwise registration. Our method is built upon the Expectation-Maximization algorithm with regularized covariance matrices and includes spatial regularization. The number of appearance regions is determined by a novel model order selection criterion. The prior is described on a reference mesh where each vertex has a probability to belong to several intensity profile classes.

Introduction

Intensity profiles were among the first image representations used to describe appearance for segmentation purposes. Cootes used intensity profiles to build Statistical Appearance Models [START_REF] Cootes | Using grey-level models to improve active shape model search[END_REF]. They are sampled in training images and both mean profile and its principal modes of variation are extracted for each landmark. Intensity and gradient profiles were used to optimize image forces of deformable models [START_REF] Gilles | Anatomical and kinematical modelling of the musculoskeletal system from mri[END_REF]. The idea was to better discriminate organ contours in images by comparing intensity profiles, using two generic models and checking the similarity variation with the normalized cross-correlation. A thorough study on intensity profiles can be found in [START_REF] Heimann | Statistical Shape Models for 3D Medical Image Segmentation[END_REF].

Several issues may be raised with statistical appearance methods based on Principal Component Analysis (PCA). First, they require an accurate point-wise registration as the statistical analysis of appearance is performed at each point. Defining homologous points for 3D structures is difficult and therefore registering those points accurately is still considered challenging. A second limitation common to most appearance models (e.g. Active Appearance Models) is that they are monomodal, i.e. they rely on the hypothesis that the probability density function is well described by a single Gaussian distribution. This hypothesis is often violated by the presence of pathologies but also by the fact that shape is not necessarily correlated with appearance (see for instance livers of Fig. 2 where top regions corresponding to the lungs vary in size). Instead of having one mode with large covariance, it is preferable for image segmentation or image detection purposes to have several modes with lower covariance.

In this paper, we describe how to build a Multimodal Prior Appearance Model from a training set of P meshes. Intensity profile classes are estimated for each mesh and not for each point (i.e. without the need for any registration). Registration between subjects is only used to estimate the posterior probabilities on a reference mesh. Furthermore, the proposed method is fully automated and one single threshold J controls the number of classes. Finally, we introduce new regularization strategies of covariance matrices in the Expectation-Maximization algorithm (EM) and the OSI index to determine the optimal number of classes.

Building Multimodal Prior Appearance Models

Principles

In Fig. 1, we overview our automated method to create a Multimodal Prior Appearance Model. The input is a set of P meshes M p corresponding to the segmentation of the same structure in different images. The meshes may have different number of vertices, or even different topologies. At each vertex i of M p , we extract M regularly sampled intensities to build an intensity profile of dimension M along the normal direction, noted x p i . This profile can extend inward, outward or both sides, depending on the application. Note that intensity profiles act as feature vectors that could be replaced by any other local or global features such as isophote curvature, texture descriptors, oriented filters, etc. Changing the feature vector would only change the regularization of the covariance matrices (section 2.2). For each mesh, we propose to automatically cluster the profiles using an EM classification. The number of classes, a hyperparameter, is selected in an automatic fashion through a model order selection based on a new criterion (section 2.3). Classification is improved by performing spatial regularization of the posterior probabilities (section 2.4). The creation of the prior is done in two steps. First, all intensity profile classes from the P subjects are compared and classes corresponding to the same tissues are possibly merged. Finally, all P meshes are registered to the same reference mesh M ⋆ and a unique prior model is created. Each vertex i of M ⋆ is given a probability γm i to belong to a reference class m (section 2.5).

EM classification of intensity profiles

The classification for a given mesh M p is formulated in the context of a probability density estimation using Gaussian Mixture Models. The approach is semiparametric because the number of classes K p is an unknown parameter. In the remainder, intensity profile class is denoted by mode and we write µ p k , Σ p k , π p k and p γ k i for the mean, covariance matrix, mixing coefficient and posterior probabilities of class k from mesh M p .

Initialization EM is initialized with the Fuzzy C-Means algorithm (FCM). After convergence, FCM cluster centers and data membership values are used to initialize the EM mode means µ p k , posterior probabilities p γ k i , mixing coefficients π p k and covariance matrices Σ p k . FCM are themselves initialized with random cluster centers.

Coping with missing data Since segmented structures may be close to the image boundaries, profiles are likely to be incomplete. In order to use those incomplete profiles in the classification, a valid dimension range M i ≤ M is determined for each profile i. In the E-step, the term exp((

x p i ) T (Σ p k ) -1 (x p i )
) is only computed for the valid dimension range and multiplied by M/M i . In the M-step, the mode means and covariance matrices are normalized by the total number of valid values for each index u of the profile (1 ≤ u ≤ M ). FCM have also been extended in similar fashion to cope with incomplete data.

Covariance matrix regularization

The EM log-likelihood maximization may lead to local maxima or degeneracy. For profiles of dimension M (typically M ≈ 10), degeneracy of covariance matrices Σ p k may occur due to the coarse sampling (typically N ≈ 4000) of this high dimensional space. We propose 3 distinct methods to regularize Σ p k based on a regularization parameter h (0 ≤ h ≤ 1). The first one is based on Spectral Regularization. The covariance matrix is diagonalized Σ p k = PΛP T and the h lowest eigenvalues are set to 1% of the highest eigenvalue, thus leading to a new diagonal matrix Λh . The inverse is then computed as h (Σ p k ) -1 = P Λ-1 h P T . This can be seen as performing PCA and filtering the matrix by discarding high frequencies. In a second approach, Diagonal Regularization, the covariance matrix is regularized towards a diagonal matrix controlled by parameter h. The u, v element of matrix h Σ p k :

h (Σ p k ) u,v = (1 -h)(Σ p k ) u,v + h δ u,v (Σ p k ) u,v (δ u,v
is the Kronecker symbol). The higher h, the more diagonal dominant the covariance matrix. This approach has been used in climate modeling to cope with missing values [START_REF] Schneider | Analysis of incomplete climate data[END_REF]. The last approach, Constant Regularization, regularizes the covariance matrix towards a constant matrix Id tr(Σ p k )/M controlled by parameter h :

h (Σ p k ) = (1-h) (Σ p k )+h Id tr(Σ p k )/M .
With a high value of h, the covariance matrix converges towards a diagonal matrix with the same variance. The choice of a covariance matrix regularization method depends on the nature of the data. We have tested the three methods on profiles and we found that Diagonal or Constant Regularization method with h = 0.9 leads to the most intuitive classification results.

Model Order Selection

The objective of model order selection is to find the number of modes that best represents the data without any under or overfitting. To this end, we propose to estimate a criterion measuring the quality of EM classification for a given number of modes K. We then keep the number of modes K opt that maximizes (or minimizes) that criterion. Several criteria based on information theory have been proposed in the literature such as the Akaike Information Criterion (AIC) and the second-order AIC (AIC c ). The often preferred EM criterion is the Bayesian Information Criterion (BIC) [START_REF] Schwarz | Estimating the dimension of a model[END_REF]. We also investigated FCM criteria such as cluster validity indices [START_REF] Kim | On cluster validity index for estimation of the optimal number of fuzzy clusters[END_REF] and the Fuzzy Vector Quantization [START_REF] Saha | A new cluster validity index based on fuzzy granulation-degranulation criterion[END_REF].

In addition, we propose a new non parametric model order selection criterion called Overlap Separation Index (OSI) inspired by Kim et al. [START_REF] Kim | On cluster validity index for estimation of the optimal number of fuzzy clusters[END_REF]. This criterion is solely based on the posterior probabilities p γ k i and penalizes the overlap between modes while encouraging their separation. More precisely, the criterion is computed as OSI = C1

C2 . The first term C 1 sums the amount of overlap 2 p γ s i /( p γ r i + p γ s i ) between the best two modes r and s for profile i (i.e. modes with the highest p γ k i ). The second term C 2 is the minimum separation between any pair of modes. The separation between pair of modes r and s is computed as the sum of 2 p γ s i /( p γ r i + p γ s i ) for all profile i being classified to mode r and 2 p γ r i /( p γ r i + p γ s i ) for all profile i being classified to mode s. Table 1 shows the performance of the different model order selection criteria with varying regularization methods (spectral method has been discarded here). The number of modes being tested varies between 2 and 10, which takes around 15 minutes for a mesh with few thousand points. For outward profiles, the expected number of modes is at least 3 (air, bones, soft tissue) while for inward profiles 2 modes are expected (parenchyma and non-parenchyma). Based on Table 1 and further analysis on 6 other liver meshes, we found that OSI criterion gave the most consistent results with a limited sensitivity to regularization methods and h parameter.

Spatial Regularization

EM does not take into account the neighborhood information of profiles. This leads to non smooth posterior probability maps p γ k i , which impairs the fusion of appearance regions. To account for the connectivity between profiles, we use the Neighborhood EM algorithm (NEM) [START_REF] Ambroise | Clustering of spatial data by the em algorithm[END_REF] since it nicely extends EM and leads to efficient computation (compared to Markov Random Field). NEM is an alternate optimization of the L functional :

L( p γ k i , π p k , µ p k , Σ p k ) = L( p γ k i , π p k , µ p k , Σ p k ) + β K k=1 N i=1 N j=1 p γ k i p γ k j v ij (1)
L( p γ k i , π p k , µ p k , Σ p k ) = K k=1 N i=1 p γ k i log(π p k G(x p i |µ p k , Σ p k )) - K k=1 N i=1 p γ k i log( p γ k i ) where G(x p i |µ p k , Σ p k ))
is the Gaussian probability density function. The former term L leads to the classical EM [START_REF] Hathaway | Another interpretation of the em algorithm for mixture distributions[END_REF] while the latter is a spatial regularization term controlled by β. The neighborhood parameter v ij sets the amount of smoothing and is non-zero only if profile i is neighbor to profile j. L functional is minimized with an alternate optimization leading to a modified E-step where the posterior probabilities are iteratively estimated until a fixed point value is reached [START_REF] Ambroise | Clustering of spatial data by the em algorithm[END_REF]. In our setup, profiles are extracted from 2-simplex meshes for which each vertex has only 3 neighbors [START_REF] Montagnat | Globally constrained deformable models for 3D object reconstruction[END_REF]. Thus, v ij has only three non-zero values, which substantially speeds-up the computation. In practice, less than 5 iterations are necessary to obtain stable posterior probabilities.

As neighborhood parameter v ij , we choose the correlation coefficient between neighboring profiles i and j. With this choice, the spatial regularization of posterior probabilities is stronger between similar neighboring profiles, similarly to anisotropic diffusion in image processing. This prevents the blurring of tissue modes that would have occurred with a constant v ij value.

The choice of the β parameter is an important issue and we set this parameter automatically by using an heuristic proposed by Dang [START_REF] Dang | Spatial fuzzy clustering using em and markov random fields[END_REF]. It consists in using NEM with increasing values of β, and then detecting the β value above which the log-likelihood L sharply decreases (see left of Fig. 3). Indeed, too much spatial regularization leads to a significantly worse profile classification captured by the L functional. The proposed approach is a fully automatic way to spatially regularize posterior probabilities.

Fusion of modes from P meshes

The objective of this section is to merge the modes from the profile classification performed on P meshes into a single Multimodal Prior Appearance Model. This is done in two steps. The first step consists in comparing and merging the appearance regions extracted on the same structure for P different subjects. Profiles of each mesh M p have been classified and lead to K p modes (K p may vary among meshes, e.g. due to the occurrence of pathologies). In order to have a meaningful comparison, an intensity normalization is required (e.g. to cope with the different nature or settings of the imaging systems). This may be done by histogram normalization or many other approaches proposed in the literature. With CT images, we found best not to perform any normalization.

In order to merge similar regions, we measure the similarity between any pair of modes (µ p k , Σ p k ), (µ q l , Σ q l ) for p = q by using the Jaccard index (ratio of the intersection of two sets over their union) of the region spanned by the mean and standard deviation µ p k ± σ p k where σ p k is the diagonal of the covariance matrix (see Fig. 4). A threshold J between 0 and 1 is used to decide whether two modes k and l are equivalent. Thus, we create a graph where nodes represent the modes and arcs link the modes found to be equivalent. The number of connected components of this graph is the number of independent modes K. For connected components with only one node (i.e. without equivalence), modes are directly included in the prior with a new index m. For connected components having r equivalent nodes, we compute the mean of the new mode as the weighted sum of profile means µ p r with the weight p i γ r i while covariance matrices are recomputed. This computation leads to K independent modes (μ m , Σm ) and an equivalence table η(p, m) establishing the new index m of mode k.

An alternative to this first step could be to perform an EM classification of all profiles for all P subjects with model order selection to find the optimal number of modes. This approach would lead to a more time consuming task, which would need to be performed each time a new dataset is added. Instead, we prefer to achieve a separate clustering of each dataset followed by a merging of all modes. Another advantage of our approach is that it is not biased by the variation of mesh resolution between datasets.

The second step provides a geometric embedding for the independent modes. To this end, we register non-rigidly all P meshes M p towards the same reference shape with a coarse-to-fine deformable surface approach where each mesh is registered towards a binary image with globally-constrained deformations [START_REF] Montagnat | Globally constrained deformable models for 3D object reconstruction[END_REF]. After defining a reference mesh M ⋆ on the reference shape, each posterior probability p γ k i is resampled on M ⋆ using a closest point approach. Finally, for each vertex i of M ⋆ , we compute the posterior probability γm i by summing and normalizing the posterior probabilities associated to each mode : γη(p,m) i + = p γ m i . In practice, this approach leads to sparse probabilities where only a few modes have non-zero posterior probabilities (as opposed to performing an E-step based on the mode means and covariances).

Results

We tested our method on 7 livers segmented from CT images and 4 tibias (cropped at knee level) segmented from MR images (see Fig. 2). For both structures, outward profiles (10 samples extracted every mm) were generated from meshes with ≈ 4000 vertices. As said before, EM classification using Diagonal or Constant Regularization with h = 0.9 leads to the most intuitive results. The optimal number of modes were estimated with the OSI criterion. NEM was launched ≈ 10 times to find the optimal β (see left of Fig. 3). For the livers, an initial total number of 24 modes leads to 14 new modes after the merging of profiles with J = 0.6. With a lower J = 0.5, the number of modes goes down to 6, thus providing a simple way to taylor the complexity of the prior. For tibias, an initial total number of 11 modes leads to 6 new modes with J = 0.5. For both structures, the Multimodal Prior Appearance Model has been built from these new modes (see Fig. 5).

Conclusion

We propose in this paper a method that builds a Multimodal Prior Appearance Model from the regional clustering of intensity profiles. The main advantage of our approach is that modes are built without requiring an accurate pointwise registration. Another advantage is that a meaningful prior may be built with very few datasets (in fact one dataset suffices), which makes it well suited for a bootstrapping strategy. Furthermore, the prior is multimodal therefore able to cope with large variation of appearances including pathologies. We have also introduced the OSI index and included spatial anisotropic regularization of EM classification. Future work will focus on the use of this prior in image segmentation.

Fig. 1 .

 1 Fig. 1. Proposed pipeline for Multimodal Prior Appearance Model construction.

Fig. 2 .

 2 Fig. 2. EM classification of outward profiles performed on 4 livers and 2 tibias.

Fig. 3 .

 3 Fig. 3. L w.r.t β (left), liver before (middle) and after (right) NEM regularization.

Fig. 4 .

 4 Fig. 4. Left: µ p k (solid lines) and σ p k (dashed lines) of 2 modes (asterisk and square) after EM classification of profiles (σ p k is the diagonal of Σ p k ). Right: Similarity between pairs of modes defined as the ratio of the intersection (dark gray) over the union (light gray) of their variance surface.

Fig. 5 .

 5 Fig. 5. Multimodal Prior Appearance Model for the livers (1) and tibias (2). For both structures, plot of the new mode means (a) and the new classification (b).

Table 1 .

 1 Selection of the optimal number of EM modes for outward and inward profiles.

			Outward profiles			Inward profiles
	Regularization h OSI FVQ AIC AICc BIC OSI FVQ AIC AICc BIC
	Diagonal	0.9 4	2	3	3	2	2	2	4	4	2
	Diagonal	1.0 3	2	3	3	2	2	2	4	4	2
	Constant	0.9 4	2	5	5	2	2	2	6	6	2
	Constant	1.0 3	2	3	3	2	2	2	5	5	2

Acknowledgment This work is supported by the EU Marie Curie project 3D Anatomical Human (MRTN-CT-2006-035763).


Hans Lamecker

Xavier Pennec

Atlas to Image-with-Tumor Registration based on Demons and Deformation Inpainting

This paper presents a method for nonlinear registration of images, where there exists no one-to-one correspondence in parts of the image. Such a situation occurs for instance in the case where an atlas of normal anatomy shall be matched to pathological data, such as tumors, resections or lesions. Our idea is to use local confidence weights and to model pathological regions with zero confidence. We integrate this concept into the efficient and publicly available diffeomorphic demons registration framework. Finally, we show that this process better captures deformations in high-confidence regions than without using the proposed modification. Furthermore, it is easy to implement and runs faster than previous approaches.

Introduction

Spatial registration of a 3D atlas to patient images is a fundamental approach to enhance individual patient data with valuable information, e.g. for simulating therapy effects. Such information may be a segmentation, physiological or structural data, which is not available for a given patient, as well as other clinical parameters. Reasons for the unavailability of such patient information in a clinical setting are manifold. Factors like radiation exposure, anesthesia, invasiveness of the measurement process, etc. have to be taken into account to find the correct balance patient safety, comfort and therapy benefit.

Atlas to patient registration of 3D image data faces a variety of challenges. First of all, the construction of an atlas is not trivial. Often, such a construction involves application-specific decisions, which renders the use of an atlas for other applications difficult. These decisions influence the choice of an appropriate image-to-image registration scheme, including cost functions (e.g. single or multi-modal), class of deformations (especially when dealing with large deformations), regularization, etc. More importantly, in the usual clinical setting patient data is pathological. One particularly challenging case are topological changes between the atlas (representing normal anatomy) and the patient data, e.g. because of the presence of tumors, lesions or surgical resections. Under these circumstances, the deformation of the atlas onto the patient image cannot be modeled with a smooth deformation.

A general solution for this problem seems beyond reach, hence applicationspecific methods are devised, which incorporate additional knowledge about the type of pathology at hand. For instance, in atlas to tumor registration one approach is to incorporate tumor growth models and seeding strategies into the registration process [START_REF] Bach Cuadra | Dense deformation field estimation for atlas-based segmentation of pathological MR brain images[END_REF][START_REF] Zacharaki | Non-diffeomorphic registration of brain tumor images by simulating tissue loss and tumor growth[END_REF]. This is prone to errors in the initial seed placement.

An alternative approach is to endow a registration method with a local measure of confidence in the resulting deformation [START_REF] Stefanescu | Non-rigid atlas to subject registration with pathologies for conformal brain radiotherapy[END_REF]. Then the pathological region can be modeled with low or zero confidence. This approach is more general in the sense that it decouples the registration from the specific pathology model.

The main contribution of this work is to combine the concept of zero-confidence local weights with the highly efficient log-domain diffeomorphic demons framework [START_REF] Vercauteren | Diffeomorphic Demons: Efficient Non-parametric Image Registration[END_REF] in such a way that its efficiency is preserved. This will be accomplished through a inpainting process as an additional step in the demons scheme (see Sec. 2). In [START_REF] Stefanescu | Non-rigid atlas to subject registration with pathologies for conformal brain radiotherapy[END_REF], the regularization step of the demons is modified over the whole image, while our approach only requires a computation on the zero-confidence region. Furthermore, our method can be implemented easily as an extension to the freely available ITK implementation of the demons [START_REF] Dru | An ITK Implementation of the Symmetric Log-Domain Diffeomorphic Demons Algorithm[END_REF].

We will also show in a prototypical study that this improves the segmentation of structures around pathologies both on simulated and clinical data sets, and can thus serve as a basis for further therapy planning.

Demons Algorithm with Deformation Inpainting

In this section we briefly summarize the demons method and motivate and explain our extension. Furthermore we compare it to related approaches.

Image registration can be formulated as a minimization problem of a cost function E(F, M ; s), that measures the dissimilarity between a (fixed or reference) image F and a deformed (moving or template) image M • s, where s is a spatial deformation. Usually, such a cost function consists of two terms: an intensity dissimilarity E i and a regularization term E r . One fundamental idea to explain the demons algorithm lies in the introduction of a hidden variable [START_REF] Cachier | Iconic Feature Based Nonrigid Registration: The PASHA Algorithm[END_REF], called the correspondences c, as well as an additional cost E c measuring the deviation between correspondences and deformation:coefficients

E(c, s) = 1 σ 2 i E i (F, M, c) + 1 σ 2 c E c (c, s) + 1 σ 2 r E r (s) (1)
where

σ i = unreliability of image intensity σ c = uncertainty of correspondences σ r = irregularity of the deformation
This idea then allows to come up with an alternate minimization scheme starting with k = 0 as follows:

1. Correspondence step: given s k , solve

c k = argmin c E 1 (c) := 1 σ 2 i E i (F, M, c) + 1 σ 2 c E c (c, s k ) (2)
2. Regularization step: given c k , solve

s k+1 = argmin s E 2 (s) := 1 σ 2 c E c (c k , s) + 1 σ 2 r E r (s) (3)
This type of regularization is called diffusion-like, while when s is replaced with ∂s/∂t in (3) it is called fluid-like. For suitable choices of E i , E c and E r the minimization can be performed very efficiently, in particular for

E i = F -M • s 2 . Choosing E c = c -s 2
and E r as a weighted sum of spatial derivatives of s leads to a Tikhonov functional, which can be solved in Fourier space. One special case are separable Gaussian kernels, which can be implemented efficiently, and are hence commonly used for regularization. For more details refer to [START_REF] Cachier | Isotropic Energies, Filters and Splines for Vector Field Regularization[END_REF]. Finally, the whole program was consistently extended to deal with diffeomorphic deformations as described in [START_REF] Vercauteren | Diffeomorphic Demons: Efficient Non-parametric Image Registration[END_REF]. Instead of replacing the global uncertainty coefficients with local coefficients as proposed in [START_REF] Stefanescu | Non-rigid atlas to subject registration with pathologies for conformal brain radiotherapy[END_REF], which leads to a highly non-trivial mathematical optimization problem for the regularization step, we propose an extension to the optimization scheme of ( 2) and (3). Let us denote the pathological region Ω in the fixed image F , and assume the atlas is given by the moving image M : If the input image is given on a regular grid, the Laplace equation turns into a sparse linear system of equations, where at each voxel i inside Ω we have

v i = 1 |N i | j∈Ni v j (4)
where v i is the local deformation at voxel i, and N i is the set of neighboring voxel to i. All v i on the boundary ∂Ω of the domain are fixed. One option would be to solve this system once and apply the inverse matrix at each iteration of the demons. However, storing this matrix can become a prohibitive burden on computer memory even for small pathological regions. An alternative way is to use a fixed point technique such as iterated means, whose stationary point is a solution of the Laplace equation. This is easy to implement: given the values on ∂Ω iteratively compute (4) for all voxels in Ω until v i does not change more than a given threshold. This approach models zero confidence in the deformation within the pathological region, since no information about the image data is considered inside Ω. The inpainting step merely interpolates the deformation on the boundary of the domain into the pathological region. As before, depending on whether we apply the step to s or the update of s we get diffusion-or fluid-like inpainting. This step also generalizes to diffeomorphic deformations by performing the inpainting in the log-domain, either on the velocity or the update field [START_REF] Vercauteren | Diffeomorphic Demons: Efficient Non-parametric Image Registration[END_REF].

The influence of undersegmentation of the tumor is different to that of oversegmentation. The latter would correspond to a relaxation of the v i on the tumor boundary, reflecting a 'not certain' region (as opposed to the 'certainly not' within the tumor).

The idea of inpainting was also presented in [START_REF] Sdika | Nonrigid registration of multiple sclerosis brain images using lesion inpainting for morphometry or lesion mapping[END_REF], but as a pre-processing step before the registration. Here the pathological regions in the patient image are inpainted directly. This may introduce a bias since arbitrary new image information is added. In [START_REF] Stefanescu | Non-rigid atlas to subject registration with pathologies for conformal brain radiotherapy[END_REF], a diffusion equation is solved with a stiffness field that locally varies across the whole image. Their method was not designed for diffeomorphic deformations, and is implemented using a grid-based parallel infrastructure. Hence, in its sequential version their algorithm runs significantly slower than ours. Another approach was presented in [START_REF] Cahill | A Demons Algorithm for Image Registration with Locally Adaptive Regularization[END_REF]. It involves solving a coupled system of PDEs. Here, the generalization to diffeomorphic demons is not obvious, and no runtimes are provided. More importantly, it is not clear how to model zero-confidence in this approach.

Results

Simulated MR Data

As a first validation we compare the result of the conventional registration (steps 1 and 2 only) [START_REF] Dru | An ITK Implementation of the Symmetric Log-Domain Diffeomorphic Demons Algorithm[END_REF] with the modified method (including step 3) on simulated T1 weighted MR data. The system presented in [START_REF] Prastawa | Simulation of brain tumors in MR images for evaluation of segmentation efficacy[END_REF] allows to place tumor seeds in a normal reference image (in this case a data set from the BrainWeb database [START_REF] Kwan | MRI simulation-based evaluation of imageprocessing and classification methods[END_REF]), and simulate their growth. The result is a simulated MR image with a tumor. However, the exact growth process is known and hence a segmentation of all structures in the tumor image, including white matter, is available.

We register the normal reference image from the BrainWeb database with the simulated tumor image (see Fig. 1a) whose segmentation around the tumor is known (Fig. 1b). Fig. 1c shows the deformation of the segmentation from the atlas image to the tumor image using the conventional registration, where there is a clear bias due to the distinct intensity in the tumor region. This bias largely disappears (Fig. 1d) when the inpainting step is switched on. Almost everywhere on the boundary of the tumor the white matter deformation is captured much more accurately, i.e. the segmentation in Fig. 1d is much more similar to Fig. 1b than Fig. 1c is to Fig. 1b. In both cases we use the same parameter settings for the demons algorithm. In particular we adopt diffusion-like regularization only. Experiments with fluid-like regularization showed inferior results.

Clinical Data

In this experiment we perform an evaluation on clinical T2 weighted images of a patient with a glioma. We use two image acquisitions of the patient at two different time-points about 5 months apart. During this time the glioma has evolved quite significantly (Fig. 2). For both images, a segmentation of the tumor is assumed to be given (in our case it was segmented in the T2 images by an expert). The atlas is a segmented T2 image of a normal subject (Fig. 3). In order to compare the registration performance with and without inpainting, we perform the following experiment. Let us denote the two patient acquisitions with M and N , and the atlas image with A. First, we register M with N to obtain a deformation T : M → N . This is a conventional (smooth) reference registration [START_REF] Dru | An ITK Implementation of the Symmetric Log-Domain Diffeomorphic Demons Algorithm[END_REF], since both images contain a tumor. Therefore we do not expect topological changes (no metastases). Next, we register the atlas A to both timepoints to obtain T 1 : A → M , and T 2 : A → N . The situation is depicted in the diagram in Fig. 4. Taking the reference deformation T between M and N as ground-truth we expect to get a similar deformation by taking the detour via the atlas. Therefore we measure the norm of these difference as

T (p) -(T 2 • T -1 1 )(p)
for each voxel p and both cases without and with inpainting. The result is illustrated in one axial slice through the image in Fig. 5. Blue-like colors denote small deviations while red-like colors indicate deviations above 5 mm. We note that the deviations around the tumor boundaries reduce significantly when the inpainting step is applied as compared to the conventional registration process. This means, that there is less bias from the image intensity inside the tumor regions in the case of inpainting deformations, hence the registration becomes more consistent with the reference deformation in the critical regions along the tumor boundaries.

Performance

The method was implemented as an extension to the ITK diffeomorphic demons code [START_REF] Dru | An ITK Implementation of the Symmetric Log-Domain Diffeomorphic Demons Algorithm[END_REF], using iterated means to solve the inpainting problem. The clinical data sets have a resolution of 256 × 256 × 64 voxels. On an Intel Xeon processor with 3 GHz the runtime for the conventional registration was about 1 min, and for the modified demons method about 9 min, where the additional time stems from the inpainting step. This time may be reduced using the inverse matrix approach on the cost of memory consumption (the glioma region consists of about 50k voxels, which would result in a matrix size of about 10 GB in single precision). For the iterated means approach, the runtime depends on the number of voxels in the pathological region and the stopping criterion for the local averaging. In our experiments, we stop when the change of the local average reaches machine accuracy everywhere. This condition may be softened to increase speed.

Conclusions and Outlook

In this paper we have proposed a simple extension of the diffeomorphic demons framework with the aim of improving atlas to patient registration in the presence of pathologies in the patient image. This extension is easy to implement, as it does not modify the demons framework but constitutes a post-processing step, which can be implemented easily. This additional step can be seen as a inpainting process where the current deformation field is interpolated from the boundary of a pathological region into the region itself. The rationale behind this approach is that the correspondences and hence the deduced deformation field is not reliable in the pathological region because there is no one-to-one correspondence with a region in the atlas.

The use of the Laplacian with fixed boundary values does guarantee a linear transition from one end of the gap to the other end. However, it does not guarantee smoothness at the boundary, meaning that the resulting inpainted vector field probably isn't diffeomorphic. One avenue for future research is to find an inpainting method that will preserve the diffeomorphic properties of the original demons algorithm. Furthermore, from a theoretical perspective it is interesting to derive the full mathematical formulation of locally adaptive confidence in the demons framework. As a next step, the connection may be analyzed between this formulation and the inpainting process as a kind of limit as the confidence tends to zero.

The inpainting procedure allows for the integration of more general (geometriclike) knowledge into the registration process by using other models than the Laplace model of this work. The image inpainting literature [START_REF] Shen | Inpainting and the Fundamental Problem of Image Processing[END_REF] offers a rich va-riety of ideas, e.g. for incorporating anisotropy to model structural or functional tissue properties.

Application-relevant issue are: (a) selection of the atlas, and (b) validation of the performance of the approach when transferring structural data such as DTI fiber tracks, for example. In pediatric oncology, the acquisition of DTI is prohibitive due to the acquisition time and consequential artifacts. One type of benchmark could be the evaluation of the therapy outcome based on a tumor evolution prediction.

3 .

 3 Inpainting step: Given the deformation s k+1 on the boundary ∂Ω of the pathological domain, overwrite the current deformation s k+1 inside Ω with the solution of the Laplace equation ∆s = 0.

Fig. 1 .

 1 Fig. 1. (a) shows simulated tumor image. (b) shows manual segmentation of white matter (turquoise) around tumor (red), and a ROI on the tumor boundary (orange ellipse). (c) shows atlas-based segmentation without deformation inpainting (bad segmentation in ROI), (d) with inpainting (much better segmentation in ROI).

Fig. 2 .

 2 Fig. 2. T2 weighted patient images with glioma at two different time-points.

Fig. 3 . 2 Fig. 4 .

 324 Fig. 3. Atlas: T2w MRI of normal subject

Fig. 5 .

 5 Fig. 5. Differences between reference deformation T and concatenated atlas deformations T2 • T1 without (left) and with (right) inpainting (color scale: blue 0 mm to red 5 mm).

Acknowledgments The authors thank Pierre-Yves Bondiau for providing clinical data. Hans Lamecker was funded by the EU-FP6 Project Health-e-Child (IST-2004-027749). Thanks to Ezequiel Geremia for helpful discussions.


E Pernod

M Sermesant

J Relan

H Delingette

Interactive real time simulation of cardiac radio-frequency ablation

Keywords: (according to ACM CCS): I, 6, 3 [Simulation and Modeling]: Applications-Medical interactive real time simulation

Virtual reality based therapy simulation meets a growing interest from the medical community due to its potential impact for the training of medical residents and the planning of therapies. In this paper, we describe a prototype for rehearsing radio-frequency ablation of the myocardium in the context of cardiac arrhythmia. Our main focus has been on the real-time modeling of electrophysiology which is suitable for representing simple cases of arrhythmia (ectopic focus, ventricular tachycardia). To this end, we use an anisotropic multi-front fast marching method to simulate transmembrane potential propagation in cardiac tissues. The electric propagation is coupled with a pre-recorded beating heart model. Thanks to a 3D user interface, the user can interactively measure the local extracellular potential, pace locally the myocardium or simulate the burning of cardiac tissue as done in radiofrequency ablation interventions.

To illustrate this work, we show the simulation of various arrhythmias cases built from patient specific medical images including the right and left ventricles, the fiber orientation and the location of ischemic regions.

Introduction

Virtual-reality based surgery simulation opens up new fields in medical training and planning of therapies. In this paper, we investigate the development of a simulator prototype for rehearsing radiofrequency (RF) ablation of the myocardium. This intervention which is performed through interventional radiology, consists in burning cardiac tissues responsible for arrhythmias. Arrhythmias are pathologies characterized by an abnormal heart beat or rhythm. The contraction and relaxation of the heart is indeed controlled by the depolarization and repolarization of an electric wave that propagates through the myocardium, from the sinus node, the natural pacemaker of the heart, to the ventricles. This electrical activity of the heart can be measured at the level of the thorax through the electrocardiogram but can also be measured by electrodes positioned inside the heart during RF ablation. The abnormal propagation of the electric wave inside the heart leads to dysfunction of cardiac mechanics causing a loss of efficiency when pumping the blood to the general circulation.

While there exists a large variety of cardiac arrhythmias, some of them, like atrial flutter or atrial fibrillation are fairly common with up to 6 million people concerned in Europe. Ventricular tachycardia (VT) is not as common as atrial fibrillation but is considered as a potentially life threatening pathology which may evolve into ventricular fibrillation and death. In case of severe and chronic atrial fibrillation or ventricular tachycardia, RF ablation is widely used to burn the pathological cardiac cells responsible for abnormal electrophysiological propagations. One source of difficulty for cardiologists when performing this intervention is to locate those cells based on the reading of multiple electrical signals acquired in the heart. In many cases, the therapist adopts a trial and error approach where he continuously checks if the electric propagation is altered by performed ablation. Furthermore, partly due the challenge in performing the intervention, a large proportion (around 50 %) of patients being ablated for atrial fibrillation experiments have a recurrence of the pathology at least 6 months after the intervention. c The Eurographics Association 2010.

Context and Objectives

To improve the success of RF ablation procedures, several research groups are developing computational models of the heart [BUM04] with the objective of predicting the optimal location of the ablation lines given some patient specific electrophysiological measurements and medical images of the heart function. In this paper, we propose a complementary objective which is to develop a simulator to train cardiologists in performing RF ablation. To be complete such simulators should also reproduce the catheter navigation performed under Xray guidance which contributes to the difficulty in performing the intervention. Several research groups [DMD * 08], [START_REF] Hughes | Real-time seldinger technique simulation in complex vascular models[END_REF] and commercial companies [Men], [Sim] are developing such endovascular navigation simulators, but in this paper, we have chosen to focus on the simulation of electrophysiology propagation which is suitable for real-time interaction. Fast and realistic simulation of cardiac electrophysiology is indeed fairly challenging since the physiological phenomenon occurring at the level of the cardiac cells involves different ion channels with various dynamics. Most electrophysiological models proposed in the literature [LR91] are way too complex to be implemented in realtime. Based on previous work, we have developed a multifront Eikonal model which captures the main characteristics of electrophysiology propagation (anisotropy, refractory period, restitution curve) at a low computational cost. Furthermore, we show on several examples that such model can be combined with a 3D user interface in order to simulate the pacing of the myocardium and the ablation of pathological cells.

SOFA Framework design

The proposed work has been implemented in an open source framework, SOFA. In this section, we overview the overall design of the framework but more details can be found in [ACF * 07]. This framework is targeting real-time modeling of deformable structures and mechanical behaviours, particularly for medical simulation.

SOFA is an open-source library mostly written in C++ and using XML script language. It can be used either as an external library, or used with its own user interface and 3D OpenGL views. The architecture main goals are to keep an high modularity implying, on one side, a lean core containing the overall architecture and prototype functions and on the other side, a large set of modules derived from the core. Thus, it allows to represent different aspects of the same object, e.g. mechanical behaviour, collision model, topology description or visual model, etc. Although all these models are connected together and synchronized, each one can be parametrized independently.

To create a data structure and link the different components (C++ module of the framework) together, scene graphs (simulation tree in XML format) are used to describe a simulation. In terms of simulation tree, nodes are used to create hierarchical levels of modeling (e.g. behaviour model, collision model, visual model, etc.), whereas a large collection of components corresponds to leaves under the nodes and covers a large set of fields (e.g. mechanical behaviour, force fields, topology, mass, etc. ). This description using a XML file allows an high flexibility in designing a simulation scene. Indeed, it is easy to replace one simple component by another. For example for the computation of forces acting on a mechanical system, one can replace a spring-mass force with a finite element force. This versatility also applies to the choice of the time integration schemes or collision detection pipelines. Changing a component does not affect other components in the scene due to the library architecture. Figure 1 shows one simple example of a simulation tree describing two simple objects, a cube and a plane. Both objects are described in terms of geometry and topology (the relationship between vertices) and have a specific representation for the collision detection. It is important to note that the XML scene graph describes implicitly the simulation tasks and their scheduling similar to a script pipeline. Then, the simulation tree is constantly traversed by visitors which call appropriate functions of the different components, carrying information and propagating them. Finally, each component has its own list of parameters which can be changed in the XML scene and at runtime.

Interactive electro-mechanical simulation

2.1. Improved Software Framework

Physical fields and their coupling

As explained before, the current framework has been mainly designed for the simulation of mechanical object interaction but not, for example, for the simulation of electrophysiology propagation and its coupling with mechanical objects. We have therefore extended the SOFA framework in order to simulate the evolution of a physical field, this field being any scalar, vectorial or tensorial quantity defined on a mesh. For instance, one may simulate the evolution of the electric potential (a scalar) on a mesh through a reaction-diffusion equation as shown in Figure 2.

A further extension of the framework has been to allow the coupling between mechanical and physical model evolution. To avoid confusion between models and to have an asynchronous behaviour, tag parameters defined as string identifiers were added to each component. Thus, all components having the tag "A" are processed by the visitor apart from components having the tag "B", whereas components without tags are shared with all visitors. As a result, we can design simulation scenes with a coupling between mechanical and electrophysiology variables. For example, Figure 2 shows a scene graph where the electric potential evolution is simulated on a deformable surface, the deformation of the surface having an influence on the diffusive process. Thus, on one hand, all components related to electrophysiology are tagged as "electro", and on the other hand, mechanical components are tagged as "meca". Note that the component MultiTag-MasterSolver is used to handle the sequential launch of visitors. Finally, some components such as the force fields components may require the knowledge of both mechanical and physical variables and a specific mechanism has been added to do so.

Scene Description

For the simulation of cardiac electrophysiology, it is important to define volumetric meshes of the heart. In addition, it is required to define physiological structures attached to those meshes. For instance, 2D and 3D regions on the myocardium are useful to identify the endocardium, epicardium... Fiber orientations at each tetrahedron are also needed. These information are directly generated from a medical image processing toolbox which combines image segmentation, surface extraction and volumetric mesh generation [TMD * 08].

Figure 3: Example of data stored on topological elements: fiber directions are defined for each tetrahedron in a heart mesh which has been cut and deformed.

In the design of SOFA, all these informations stored in mesh files are loaded by a MeshLoader component. This last one will store all the data into specific C++ containers, called "Data", which are C++ components members reachable in the XML file. Thus, mechanisms are available to create dependency between those containers of same nature in the XML file (similar to "fields" in OpenInventor). Moreover, some containers have been specialized to store information associated with topological mesh elements like vertices, triangles, tetrahedra... For instance, TriangleData are specific Data containers storing information on triangles (such as colour, conductivity...). This approach allows to provide an automatic mechanism to handle topological changes [START_REF] André | Versatile design of changing mesh topologies for surgery simulation[END_REF]. When simulating the ablation of cardiac tissues, the c The Eurographics Association 2010. content of those containers is automatically updated by callback functions. As example, in Figure 3, fiber orientations are updated during the resection of the myocardium.

Electrophysiology model: Fast-Marching method

Modeling the cell electrophysiology is an active research area since the seminal work of Hodgkin and Huxley [HH52]. We can divide these models into three categories, from the more complex to the simpler (numerically): Solutions of the evolution PDE are very time-consuming, due to the space scale of the electrical propagation front being much smaller than the size of the ventricles. The motion of the front governed by the Eikonal equation is observed on a much larger scale, resulting in faster computations.

To simulate arrhythmias and radiofrequency ablation, we use a multi-front Eikonal model which was first proposed in the context of electrophysiology parameter estimation [SKD * 07]. This model makes the hypothesis that at each node of the mesh, the surrounding cardiac tissue may be in one of the three successive states : depolarized, refractory and repolarized. When a node is in refractory state, it cannot be depolarized (excited) anymore. When it is in repolarized state, it can be depolarized. The duration between the depolarized and refractory states is called the Action Potential Duration or APD whereas the duration between the refractory and repolarized states is called the refractory period. The APD and refractory periods change over time and space but there exists simplified relations called the restitution curves that can be measured experimentally and indicates how those durations vary as a function of frequency.

The transition between depolarized and repolarized states takes into account the anisotropic propagation of an electrical wave governed by the following eikonal equation :

F √ ∇T T D∇T = 1
where F is the conduction velocity, T is the time at which the wave reaches a given point and D is a 3 × 3 symmetric tensor defining the anisotropy of propagation. The conduction velocity F is not homogeneous in the myocardium. Ischemic or infarcted (necrotic) tissues in particular show low conduction velocities compared to normal tissue. The anisotropy matrix D captures the fact that the electrical propagation is roughly three time as fast along the fiber orientation than in its orthogonal direction. Therefore, this tensor is computed as :

D = 1 9 I + 8 9 f ⊗ f
where I is the identity matrix and f is the unit fiber orientation vector.

To establish in which state a node should be, we use an anisotropic fast marching method [KSC * 07] which computes the depolarization time at nodes located next to the front location. One issue with anisotropic fast marching is that the front propagation direction is no longer along the gradient of the isochrones which makes the classical fast marching method [Set99] fail. The anisotropic fast marching has been implemented on tetrahedral meshes and satisfactory results in terms of accuracy and computation times even in the case of large anisotropy ratios have been obtained.

To summarize, the proposed electrophysiology model is both discrete for the state of cardiac tissue and continuous for the propagation of the electrical wave. However due to the anisotropic Fast Marching Method (FMM), it is possible to simulate the cardiac electrophysiology in real-time for tetrahedral meshes of coarse medium resolution.

Mechanical model

The electrophysiological signal propagating through the heart is responsible for the contraction and relaxation of the heart muscle, the myocardium, during the cardiac cycle. The electro-mechanical coupling of the heart is rather complex phenomenon, but has been modeled in a fairly simplified way. Several research groups [START_REF] Guccione | Theory of Heart: biomechanics, biophysics, and nonlinear dynamics of cardiac function[END_REF], [START_REF] O'sullivan | Real-time interactive volumetric animation of the heart's electrical cycle from automatically synchronized ecg: Research articles[END_REF] have developed beating heart simulators that are coupled with electrophysiology models. However, simulating the cardiac mechanics in real-time (1 cardiac cycle per second) is still currently out of reach.

As an alternative, we propose to simulate offline the heart mechanical deformation during a cardiac cycle and to store in a file the position of the mesh vertices at various instants during the cardiac cycle. During the electrophysiology simulation, the heart mesh is deformed by interpolating between the key frames stored in those files. The cardiac mechanics is simulated using the same tetrahedral mesh as for the electrophysiology modeling. This electromechanical model [SPC * 08] of the heart relies on linear tetrahedral finite elements and anisotropic linear elasticity and is able to capture the four phases of the cardiac cycle : contraction, isovolumetric contraction, relaxation and isovolumetric relaxation.

Simulation of radiofrequency ablation in case of arrhythmia

Basic User Interaction

As explained previously, FMM is an efficient computational method to approximate the propagation of transmembrane c The Eurographics Association 2010.

potential in the heart and to visualize that potential as a coloured texture. Through the topological Data containers described in section 2.1, the user can easily read and change parameters that are stored for each vertex, edge or tetrahedron. For instance, the physiological states, action potential duration are stored for each vertex while the fiber orientation and anisotropy ratio are stored for each tetrahedron. Also, the mouse can be used to change locally the conductivity of a region of the myocardium or perform an electrical stimulation to generate additional electrophysiological waves.

3D User Interface

To improve the realism of the simulation, we have interfaced our electrophysiology simulation with a 3D User Interface.

Ideally, the user should interact with a force-feedback interventional radiology simulator which would allow the user to navigate a catheter through blood vessels and perform various actions (measuring potentials and ablating tissue). However, in a first prototype, we are using a Xitact ITP instrument [Men] as 3D positioning device which consists of one tracked instrument with a button and two pedals. In fact, the ITP instrument is designed for laparoscopic simulation with force-feedback. As the appropriate tools should has been the Xitact VIST, this choice of materiel lay on logistical choice and the fact that this tool allows user interactions. For this instrument the position and orientation of its tips can be read through a serial port at a high frequency. A complete description of the handling of haptic device, positioning and force feedback in the SOFA framework, can be found in [SDC08].

Instrument Modeling

Due to nature of the 3D user interface, we assume that catheters are rigid instrument and avoid simulating the flexibility of a catheter that would be used for this kind of intervention. Figure 4 shows that the instrument is modeled as a thin rigid shaft with a sphere in its extremity. This tip is used to simulate the catheter interaction (potential stimulation and ablation). The sphere serves as a collision model for the collision detection with other objects.

In terms of XML scene graph, the instrument is a standalone node with several branches corresponding to different types of models. It consists of two mechanical and visual models, for the shaft and the tip, linked together with an articulation. It also includes a sphere as a collision model and finally a device plugin which is in charge of communicating with the hardware system. Thus, at each time step, the plugin is in charge of refreshing the 3D position and orientation of the tool. This information is then propagated through the scene graph in order to update the different components.

For further work, as indicated in previous part, the Xitact VIST model should be used for this simulation. In termes of integration in our simulator, the main part will be to simulate the catheter navigation taking into account blood fluid, vessel walls interaction and X-ray view which correspond to reality. However, some works on catether navigation has already been done in the SOFA framework [DMD * 08]. Then, the interactions developed for the ITP model should be easily transfert to this one.

Collision Detection

A sphere collision model is used to represent the tip of the tool whereas tetrahedra are employed as the collision model for the myocardium, as shown in Figure 5. Thus, as soon as the sphere intersects a tetrahedron face, a collision is created and sent to the default collision pipeline of the SOFA framework [ACF * 07]. Upon the detection of a collision, no force-feedback is sent back to the device in such a case but instead we initiate some actions through a component called the HeartSimulation-Manager. When activated by a collision, this component allows the user to perform ablation or electrical stimulation by pushing on the tool button. Also, the pedal is used to change the type of interaction. To guide the user, the selected area after the collision is highlighted and its scale may be modified as seen in Figure 5.

Results

Simulation Scene

In order to illustrate the radiofrequency ablation simulation, we build a simulation scene corresponding to ischemic ventricular tachycardia, an example of cardiac arrhythmia. In this pathology, the myocardium includes several necrotic regions or scars that appear after an infarct. Those scars can be imaged with a specific MR protocol, late enhancement MR, and are known to correspond to regions of very low conduction velocity. In the scene, we have decided to create synthetically an isthmus between two scars regions, which occurs frequently on patients suffering from VT. In this isthmus, represented in red in Figure 6, the propagation of electrical waves is also slowed down. To prevent fast heart beats and reentries, radio-frequency ablation has been shown to be an efficient therapy if performed properly [HUKS07]. During this procedure, catheters are used to permanently damage cardiac tissues with electromagnetic radiations having low voltage and frequencies in the radiofrequency range. When heated, the damaged tissue is no longer able to conduct electrical waves preventing the arrhythmia from occurring.

To simulate this pathology, a volumetric mesh of the right and left ventricles with scars has been generated from MRI images. This mesh is approximately composed of 10k vertices and 50k tetrahedra. Figure 6 shows the conductivity map on this mesh. Healthy tissues have a conduction velocity of 50 cm.s -1 whereas the scar regions (in blue) have no conduction velocities and the isthmus (in red) has a conduction velocity of 30 cm.s -1 .

Simulation Results

To display the electrophysiology simulations, we use a colour map from blue to red, to represent the electrophysiological state of each vertex. Depolarized vertices are set to red, repolarized vertices to blue and refractory vertices to yellow. At the initialisation stage, all vertices are set to the repolarized state. Then, a stimulation is performed periodically at the apex of the right ventricle at a period of 300 ms to simulate a procedure called VT Stim. This procedure is performed to test the inducibility of ventricular tachycardia of a patient and decide whether a defibrillator should be implanted. The electrical wave resulting from the stimulation takes approximately 250 ms to traverse the two ventricles. These simulations are performed with a time step of 1 ms and we obtain a mean frame rate of 45 FPS.

Figure 7 shows the resulting wave propagation. As expected, the scars remain blue, since it is not conductive whereas the isthmus zone has a different behaviour than neighbouring healthy tissues. Indeed, the wave enters the isthmus through the upper side (point A in figure 7) and then exits by its lower bound (point B). Then, due to different parameter values of the restitution curve, these vertices have a longer refractory period than healthy tissue. As a result, when the second wave (wave T+1) reaches the upper bound of the isthmus, it can no longer go through this path. However, as the isthmus is no longer in refractory period when the second wave reaches the point B, it can then enter the isthmus through this point. This pattern is called a unidirectional block, and it can generate a reentry wave.

In ischemic VT, this reentry wave can turn around the scars through the isthmus and then create additional ventricular contractions leading to tachycardia. Non-sustained VT patterns (single reentries) could be simulated with specific patterns of stimulations. In those cases, when the wave T+1 goes through the isthmus from the point B, on the other side (point A), healthy tissues are already in repolarized state. Thus, when reaching the point A, the wave will be able to propagate into the heart, as it is the case in the right part of Figure 8, and potentially sustain itself by repeating this pattern.

A second set of tests has been performed to evaluate the c The Eurographics Association 2010. importance of the anisotropy in those wave propagations. To this end, we modified the anisotropy ratio of conduction velocity in the anisotropy matrix D. To quantify those changes, we used percentages, 0% meaning that we are in an isotropic case and 90% for realistic behaviour where the electrical propagation is roughly three times as fast along the fiber orientation than in its orthogonal direction. As expected, different anisotropy ratios lead to different conduction patterns. For instance, the left part of Figure 8 shows an isotropic potential propagation whereas the right figure shows a propagation with 90% anisotropy. In the anisotropic case, the propagation includes a reentry pattern while in the isotropic case no such pattern is observed. The anisotropic Fast Marching algorithm achieves a good accuracy even in strongly anisotropic cases.

Simulating radio-frequency ablation

Additionally to the simulation described above, we simulate interactively the radio-frequency ablation therapy. In this casd, the aim is to close the isthmus in order to prevent reentry loops. Figure 9 is a screenshot of the interactive simulation while positioning the haptic device on the isthmus zone. The highlighted zone corresponds to the region impacted by the instrument, with green tetrahedra corresponding to healthy or isthmus tissues and red tetrahedra to scar tissues (with null conductivity). Adding collision detection and the prerecorded heart motion requires some additional computational time. Thus, we obtain for this simulation with a time step of 1 ms a mean frame rate of 15 FPS which is a bit lower than previous results but is still acceptable for interactive simulation.

With this setting, several radio-frequency ablations are performed around the isthmus in order to burn its tissue and get a large region with zero conduction velocity, as shown in blue in Figure 10, which blocks electric waves. To ensure that the path is completely blocked, stimulation can be performed during the intervention by actioning a pedal. By applying stimulations near the scars, one creates a new potential wave like in Figure 11 to check that the isthmus has been successfully ablated. A video clip describing the different stages of the simulation is provided as an electronic appendix.

Discussion

To assess the accuracy of the simulation, it is possible to compare the simulated isochrone map with the one measured by catheters during the intervention using contact (i.e. the NavX or Carto systems) or non contact (such as the Ensite balloon) endocardial mappings. Further evaluation can be performed by comparing stimuli conditions for which one can induce VT. In term of haptic simulation, no validation from clinicians can be performed at the moment since our prototype does not reproduce virtual catheter guidance to the heart, as it uses laparoscopic instruments for force feedback.

Conclusion and perspectives

In this paper, an electrophysiology heart model based on an anisotropic multi-front Fast Marching Method has been presented. It provides a suitable means to represent in real-time the electrophysiological cardiac rhythm and its pathologies. This model has been coupled with pre-recorded realistic heart mechanical deformations and has been integrated in the open source framework SOFA, for which specific software components have been developed to handle the coupling between mechanical and electrophysiological behaviour.

Using a 3D user interface, we have presented a complete pipeline for the interactive simulation of cardiac radiofrequency ablation. A first prototype has been demonstrated for ischemic ventricular tachycardia where reentries may occur through an isthmus surrounded by scar regions. A 3D rigid instrument allows the user to interactively simulate radiofrequency ablation by zeroing the conduction velocity around intersected tetrahedra. The user can also apply focal stimulations to check the effect of ablations.

There are several ways to improve the realism of the current prototype. For instance, the electrophysiology model could be extended to better represent the pathological behaviour and to include patient specific restitution curves. Real-time coupling between mechanical deformations and electrophysiology propagations should be developed in order to provide additional realism. Finally, catheter navigation should be made more realistic by using interventional radiology force-feedback simulators.

Figure 1 :

 1 Figure 1: Scene graph and view of a simple system: a cube on a floor. They both have mechanical behaviours and collision models.

Figure 2 :

 2 Figure 2: Simplified scene graph and view of a coupled model: a soft tissue model with electrical potential diffusion. Each behaviour is describe by one node and a shared node where components without tags are used by both behaviour.

Figure 4 :

 4 Figure 4: Picture of the ITP Xitact tool and its corresponding model. The tool tip corresponds to a sphere collision model.

Figure 5 :

 5 Figure 5: Collision model representations. A sphere model is used for the tool whereas the myocardium collision model consists of tetrahedra. The contact region is shown in green.

Figure 6 :

 6 Figure 6: Conduction velocities of the myocardium. Zones in wireframe have normal conductivity while blue regions (scars) have very low values and red region (isthmus) has intermediate values.

Figure 7 :

 7 Figure 7: Electrophysiological wave propagation simulated with a multi-front Eikonal model. Successive waves are represented (T and T+1) as well as isthmus boundaries (points A and B).

Figure 8 :

 8 Figure 8: Potential wave propagations with different anisotropy ratios; (Right) isotropic propagation; (Left) 90% anisotropy propagation with a reentry wave at the top boundary of isthmus created by wave T+1.

Figure 9 :

 9 Figure9: Simulation of radio-frequency ablation on the isthmus performed with a 3D user interface. The green zone corresponds to the tissues to be burned whereas the red zone represents tissues with no conductivity (scars).

Figure 10 :

 10 Figure 10: Wireframe representation of conduction velocity which reveals a large scar region (blue volume) resulting from the radio-frequency ablation.

Figure 11 :

 11 Figure 11: Potential waves propagation simulated after an ablation of the isthmus. Using the 3D user interface a stimulation is applied to generate a second wave near the isthmus.

c The Eurographics Association 2010.