Integration of electro-anatomical and imaging data of the left ventricle: An evaluation framework

Medical Image Analysis 32 (2016) 131–144

Contents lists available at ScienceDirect

Medical Image Analysis journal homepage: www.elsevier.com/locate/media

Integration of electro-anatomical and imaging data of the left ventricle: An evaluation framework David Soto-Iglesias a,∗, Constantine Butakoff a, David Andreu b, Juan Fernández-Armenta b, Antonio Berruezo b, Oscar Camara a a b

PhySense, DTIC, Universitat Pompeu Fabra, Barcelona, Spain Arrhythmia Section, Cardiology Department, Thorax Institute, Hospital Clinic, Barcelona, Spain

a r t i c l e

i n f o

Article history: Received 21 July 2015 Revised 22 January 2016 Accepted 25 March 2016 Available online 4 April 2016 Keywords: Quasi-conformal mapping Electro-anatomical mapping MRI Cardiac arrythmias

a b s t r a c t Integration of electrical and structural information for scar characterization in the left ventricle (LV) is a crucial step to better guide radio-frequency ablation therapies, which are usually performed in complex ventricular tachycardia (VT) cases. This integration requires finding a common representation where to map the electrical information from the electro-anatomical map (EAM) surfaces and tissue viability information from delay-enhancement magnetic resonance images (DE-MRI). However, the development of a consistent integration method is still an open problem due to the lack of a proper evaluation framework to assess its accuracy. In this paper we present both: (i) an evaluation framework to assess the accuracy of EAM and imaging integration strategies with simulated EAM data and a set of global and local measures; and (ii) a new integration methodology based on a planar disk representation where the LV surface meshes are quasi-conformally mapped (QCM) by flattening, allowing for simultaneous visualization and joint analysis of the multi-modal data. The developed evaluation framework was applied to estimate the accuracy of the QCM-based integration strategy on a benchmark dataset of 128 synthetically generated ground-truth cases presenting different scar configurations and EAM characteristics. The obtained results demonstrate a significant reduction in global overlap errors (50–100%) with respect to state-of-the-art integration techniques, also better preserving the local topology of small structures such as conduction channels in scars. Data from seventeen VT patients were also used to study the feasibility of the QCM technique in a clinical setting, consistently outperforming the alternative integration techniques in the presence of sparse and noisy clinical data. The proposed evaluation framework has allowed a rigorous comparison of different EAM and imaging data integration strategies, providing useful information to better guide clinical practice in complex cardiac interventions. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Ventricular tachycardia (VT) caused by re-entry circuits is one of the most critical forms of arrhythmia, which is usually treated with anti-arrhythmic drugs and, in more severe cases, with radiofrequency ablation (RFA) (Zipes et al., 2006). The aim of RFA is to eliminate the myocardial substrate responsible for the VT. It is generally performed using a special navigation system such as CARTO (Biosense, Cordis Webster, Marlton, NJ), which provides an electroanatomical voltage map (EAM) that characterizes the electrical be-

∗

Corresponding author. Tel.: +34699514271. E-mail addresses: [email protected] (D. Soto-Iglesias), [email protected] (C. Butakoff), [email protected] (D. Andreu), [email protected] (J. Fernández-Armenta), [email protected] (A. Berruezo), [email protected] (O. Camara). http://dx.doi.org/10.1016/j.media.2016.03.010 1361-8415/© 2016 Elsevier B.V. All rights reserved.

havior of the heart. In these contact-mapping systems a EAM is reconstructed from the processing of unipolar or bipolar electrograms acquired through the contact of a catheter with the endocardial (or epicardial) wall in a discrete set of points. Peak-to-peak amplitudes of 1D electrograms are estimated to generate both Local Activation Time (LAT) and voltage maps. Such maps are used to identify both the myocardial substrate and potential ablation targets by finding electrogram abnormalities such as low voltage regions, delayed activation times or double potential events. However, the success rate of RFA is still low (recurrence rates up to 91% (Berruezo et al., 2012) in some patients) and the ablated area is usually larger than optimal, partially due to some limitations of the EAM acquisition system: availability of a reduced number of measurements, lack of transmural information, the effect of farfield signals, uncertainties in the EAM point localization or imperfect catheter contact among other reasons.

132

D. Soto-Iglesias et al. / Medical Image Analysis 32 (2016) 131–144

In order to improve the localization of the ablation targets and to reduce the ablated areas in scar-related VTs, advanced ablation strategies have been developed (Fernández-Armenta et al., 2013). Such VTs are commonly observed in post-myocardial infarction patients, which present heterogeneous scars composed by a nonconducting zone, known as core zone (CZ), and a slow-conducting one, often referred to as border zone (BZ) or peri-infarct zone. Border zone regions are extremely important since they can be located within a non-conducting zone linking two healthy tissue (HT) regions, thus creating a conduction channel (CC) that forms a circuit of action potential re-entry. One of the most successful ablation strategies, the so-called dechannelling technique (Berruezo et al., 2012), precisely aims at only ablating the entrances of the CCs, eliminating thus the re-entry circuit responsible for the VT with a minimal damage of the BZ and HT regions. Nevertheless, identifying CCs exclusively from EAMs requires considerable operator experience, suffers from substantial inter- and intra-observer variability and is highly time consuming (Codreanu et al., 2008). Some studies have shown how delay-enhancement magnetic resonance images (DE-MRI) can assist on characterizing the scar (BZ and CZ) prior to the intervention, identifying possible locations of the CCs and incorporating the extracted information into clinical VT therapy planning (Andreu et al., 2011). The BZ channels are defined as continuous corridors of BZ surrounded by scar core/mitral annulus (Fernández-Armenta et al., 2013). A corridor of BZ is then considered a channel when connecting two areas of normal myocardium. Relating CC information obtained from EAM and DE-MRI is challenging due to the intrinsic limitations of the modalities (e.g. limited number of slices and anisotropic resolution in DE-MRI) and their different acquisition nature (e.g. coordinate systems, different and/or incomplete fields of view). A robust approach is then needed for establishing correspondences between these complementary sources of information. In clinical routine correspondences are usually obtained in two stages1 : first applying a rigid registration technique guided by a set of landmarks manually selected in both modalities; then mapping information available at each point of the EAM surface onto the closest point of the LV geometry segmented from DE-MRI. During the registration stage, point-to-surface distance and landmark alignment are both optimized to determine the best transformation (LiFern, 2008). Nevertheless, a registration technique only based on landmarks might not be sufficient to compensate for spatial differences between EAM and DE-MRI data, as demonstrated by a recent study (Andreu et al., 2015), where 79.2% of all the EAM CC’s could be visually identified in the corresponding DE-MRI data. Some researchers (Tao et al., 2012; Roujol et al., 2012) have incorporated scar information in the registration process to further constrain rigid transformations between EAM and DE-MRI data. This strategy strongly depends on the scar definition in both modalities, which in turn is highly influenced by the thresholds employed to classify scar tissue types (Andreu et al., 2011). Furthermore, all these advanced registration techniques are hampered by subsequently using a simple closest point mapping method once EAM and DE-MRI surface meshes are registered. A thorough comparison among these different integration techniques cannot be easily performed due to the lack of ground-truth data in this application since EAM and DE-MRI data represent different, even if related, physical phenomena. For instance, scar regions are identified from voltage amplitude of the EAM signals, while DE-MRI gives information about how much time the injected Gadolinium contrasts has been infiltrated in myocardial tissue (necrotic tissue taking more time to release the Gadolinium

1 Within the CARTO system (CARTO, Biosense, Cordis Webster, Marlton, NJ), correspondences are obtained with the software CartoMerge.

out). It is quite obvious that both modalities will provide different scar definitions, independently of the integration strategy. In this paper, we propose a rigorous evaluation framework, with appropriate global and local measures, to assess the accuracy of the integrated multimodal information, allowing a detailed and quantitative analysis of the process as a whole, but also to identify which steps (e.g. registration, mapping) and parameters are the most critical for the final result. Some of the developed evaluation measures are particularly tailored for EAM and DE-MRI data of the left ventricle (LV) since they relate EAM critical sites (CC entries, double potential points) with scar tissue information (CZ, BZ, HT) extracted from DE-MRI. We also propose a new mapping technique for integrating EAM and DE-MRI information that relies on a conformal mapping between LV endocardial surfaces and a 2D disk, followed by a correction of anatomical landmarks based on Thin-Plate Splines (TPS) that relaxes the established conformal mapping to quasi-conformal (QCM). Several researchers have recently proposed mapping techniques to construct 2D reference systems for different organs such as the left ventricle (De Craene et al., 2012), the atria (Karim et al., 2014; Tobon-Gomez et al., 2015), the brain (Auzias et al., 2013; Yushkevich et al., 2006), the liver (Vera et al., 2014), vertebral bones (Lam et al., 2014), the cochlea (Vera et al., 2015) or even faces (Gu et al., 2010). Information derived from both EAM and DE-MRI of a single patient at different time points or from different patients can easily be represented on the disk and then jointly be analyzed in the constructed common reference space. As an example, this approach was initially applied to compare and perform statistical analysis on different EAMs of porcine data (Soto-Iglesias et al., 2013). The developed evaluation framework was used to estimate the accuracy of the proposed QCM-based mapping technique and compare its performance to alternative state-of-the-art EAM and DEMRI integration strategies. A simple rigid transformation guided by landmarks, the Iterative Closest Point (ICP) (Besl and McKay, 1992) and a non-rigid registration technique based on currents (Vaillant and Glaunes, 2005) were analyzed for the registration stage. For the mapping stage, the commonly used closest point strategy was compared with the new quasi-conformal mapping (QCM) technique. All integration approaches were tested on a benchmark dataset of 128 synthetically generated ground-truth cases presenting different scar configurations and EAM characteristics, as well as in seventeen clinical datasets. The paper is structured as follows. In Section 2 the acquisition protocols to obtain EAM and DE-MRI data are detailed. Section 3 is devoted to the proposed evaluation framework, including the extraction of relevant information from multimodal data, the generation of the ground-truth data and the global and local measures used to estimate integration accuracy. The QCM-based mapping technique is described in Section 4, whereas the alternative stateof-the-art integration strategies are listed in Section 5. Results obtained from applying the developed evaluation framework to the different integration strategies are presented in Section 6 and finally discussion and conclusions are given in Section 7. 2. Data acquisition In this work, we applied different integration techniques to the data of 17 patients that underwent VT ablation at Hospital Clínic de Barcelona. All subjects had a DE-MRI examination prior to the VT ablation procedure using a 3T clinical scanner (Magnetom Trio, Siemens Healthcare). The 3D slab of images was acquired in the transaxial direction. Slice thickness was 1.4 mm, with no gap between slices. The field of view was set at 360 mm and matrix size was kept to 256 × 256 pixels in order to yield an isotropic spatial resolution of 1.4 × 1.4 × 1.4 mm, giving a set of images (typically

D. Soto-Iglesias et al. / Medical Image Analysis 32 (2016) 131–144

133

Fig. 1. Framework to evaluate the accuracy of EAM and DE-MRI data integration, which involves the use of clinical and synthetically generated ground-truth data. The left ventricle (LV) is segmented from DE-MRI and myocardial tissue is characterized into healthy (HT, pink), core zone (CZ, red) and border zone (BZ, green) labels, before building maps to include transmural information. The original EAM maps are re-formatted to generate a surface mesh, which is then processed for tissue characterization (HT, CZ, BZ). The synthetic data is generated from a single LV geometry from DE-MRI data, incorporating different scar configurations and noise levels. Subsequently the different integration strategies, with several alternatives for the registration (DE-MRI transmural maps and EAM wireframe mesh) and mapping stages, are applied to the clinical and synthetic data to evaluate their accuracy with global and local measures. ICP: Iterative Closest Point; QCM: Quasi-Conformal Mapping. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

50–70) in the LV short-axis orientation. In order to compensate for the expected long acquisition time, we added 80 ms to the nominal value necessary for null normal myocardium as derived from a TI-scout sequence, typically at 320 ms. Other typical sequence parameters were as follows: fat suppression, repetition time 2.6 ms, echo time 0.9 ms, flip angle 15° , read-out every other heartbeat, bandwidth 810 Hz/pixel, and 45 k-space lines filled per heartbeat in a Cartesian trajectory and antero–posterior phase-encoding direction. Image acquisition was ECG-gated to end-diastole to minimize cardiac motion. Respiratory synchronization was performed for every other heartbeat using a crossed-pair navigator approach. The dataset was acquired during expiration and parallel imaging (generalized autocalibrating partially parallel acquisition, GRAPPA) with an acceleration factor of 2 to speed up data acquisition. The mean acquisition time was 16 ± 8 min. Endocardial contact mapping data was collected with a CARTO 3 system (CARTO, Biosense, Cordis Webster, Marlton, NJ) after introducing a quadripolar catheter through the femoral vein until reaching the endocardium. EAMs were subsequently generated by reconstructing information obtained from the set of sparse acquisition points where the clinician placed the catheter to measure the electrical activity (around 2500 ms at 1 kHz). In this paper, the EAMs are reconstructed from bipolar signals, which measures the electrical potential between two electrodes of the catheter. The electrical measurements were then processed and rendered in 3D to display relevant information such Local Activation Time (LAT) or maximum bipolar intensity points (BIP) maps where heterogeneous tissue regions can be distinguished. Electrograms with delayed components (E-DCs) in a scar zone (0.5 to 1.5 mV voltage) were tagged and classified as entrance or inner CC points, depending on delayed-component precocity during sinus rhythm (Berruezo et al., 2015). The CC entrance was defined as the EDC with the shortest delay between the far-field component of healthy/BZ muscle (low frequency, usually high voltage) and local component (delayed, high frequency, usually fractionated and low voltage) corresponding to the local activation of myocardial fibers in the scar. In this study, the average number of acquired CARTO points for each patient was 555 ± 348, with a non-homogeneous

spatial distribution, more than half of these points being placed in scar-related areas. 3. Framework to evaluate the accuracy of EAM and imaging integration strategies The developed evaluation framework is illustrated in Fig. 1, showing the different steps required to pre-process the clinical data, to generate the synthetic ground-truth data and the use of global and local measures to assess the accuracy of the integration achieved by the different techniques. The first stage for processing clinical data consist in classifying the myocardial tissue into HT, BZ and CZ from both DE-MRI and EAM acquisitions, based on standard thresholds currently used in clinical routine (Andreu et al., 2011). Prior to this, a volumetric segmentation of the LV is performed in the DE-MRI followed by the generation of several 3D surface layers through the myocardium (from the endocardium to the epicardium), as in (Fernández-Armenta et al., 2013), whereupon the intensities of the DE-MRI are projected. The resulting transmural information is subsequently integrated into the endocardial layer, independently for each segmented tissue. The synthetic data generates a benchmark dataset of 128 simulated EAM maps from a single LV geometry, incorporating different scar configurations and noise levels. The different integration strategies are then applied to the clinical and synthetic data to evaluate their accuracy with global and local measures. Several alternatives for the main stages of the integration strategies are tested: (i) registration step between the DE-MRI (transmural maps) and EAM (wireframe mesh) based on landmarks, Iterative Closest Point (ICP) and currents; (ii) mapping step with the classical Closest Point (CP) approach and the Quasi-Conformal Mapping (QCM) proposed in this work, which is based on an intermediate step projecting the data into a disk representation. 3.1. Extraction of ventricular geometry and substrate information The LV geometry was extracted from DE-MRI using a segmentation algorithm based on an initial alignment of a cardiac atlas (Hoogendoorn et al., 2013) guided by manual landmarks and

134

D. Soto-Iglesias et al. / Medical Image Analysis 32 (2016) 131–144

Fig. 2. Top: three transmural layers through the myocardial wall (from endocardium, 10% layer, to mid-wall, 50% layer) with the tissue characterization provided by standard thresholds in DE-MRI data (CZ (red), BZ (green), HT (pink)). Bottom: transmural maps of tissue labels (red and pink colors represent low and high probabilities, respectively) estimated from four transmural layers (from endocardium to mid-wall). CC: conduction channel. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

followed by manual corrections where needed, as it is usually performed in clinical routine at Hospital Clínic de Barcelona. The whole myocardium was then divided into 5 transmural layers (10%, 25%, 50%, 75%, 90% of the local wall thickness starting from the endocardium), where DE-MRI pixel intensities were subsequently projected. The layer extraction was performed using a research software (ADAS, Galgo Medical SL) available at Hospital Clínic de Barcelona. In order to characterize the scar and to identify the conduction channels (CCs) of viable tissue responsible for some VTs (Andreu et al., 2011), cardiac tissue is usually classified into three types: core zone (CZ, non-conducting scarred tissue); border zone (BZ, area close to core zone with low levels of conductivity); and healthy tissue (HT). Applying this tissue classification to DE-MRI, clinicians can identify more specific potential ablation targets such as the entrances of CCs (Berruezo et al., 2015). The CZ, BZ and HT regions are identified in DE-MRI using the following thresholds on voxel intensities, as proposed in Andreu et al. (2011):

CZ ≥

60 ∗ MIV 40 ∗ MIV > BZ ≥ > HT , 100 100

(1)

where MIV is the maximum intensity value. Using nominal 60–40% thresholds, the average percentage of scar (CZ + BZ) mass with respect to the overall LV mass (tissue volume) from DE-MRI in our patient population was of 29.77 ± 11.90% (mean ± SD), which is in good agreement with similar studies such as in FernándezArmenta et al. (2013), which reported 34 ± 17.4%. A three-dimensional tissue characterization (including CC distribution) is available from the transmural DE-MRI layers but it needs to be integrated into a single (endocardial) layer to be compared with EAM data, which is only available at the endocardium. Additionally, scar transmural information can give useful insight for intervention planning since it has been associated with the location of critical sites of post-infarct VT that are potential ablation candidates (Sasaki et al., 2012). For this reason we computed three transmural maps, one for each type of tissue (HT, CZ and BZ), where each vertex of the endocardial layer is associated with a transmural probability according to the following equation:

T rans PX :

i=Endo Xi

N

(2)

where X is the type of tissue (CZ, BZ, HT) and i indicates one of the available N layers, that goes from endocardium (Endo) to the a given transmural layer (Trans). Since our database is mainly composed by subendocardial scars, we considered N = 4 layers from the endocardium to the 50% of the total wall thickness. This process then generates one layer approximately every 1.5 mm transmurally, which is in the order of the imaging resolution. The three transmural maps (HZ, CZ and BZ) provide complementary information, as shown in Fig. 2: the HT transmural map clearly separates healthy from non-healthy tissue (pink and red, respectively, in Fig. 2d); the CZ transmural map emphasizes highly transmural non-conductive scar regions (pink in Fig. 2e); the identification of conduction channels of BZ requires a combined inspection of CZ and BZ transmural maps, which can partially be visualized in both maps, i.e. high probability areas in the BZ transmural map that connect two HT regions going through a CZ area in the CZ transmural map; additionally the BZ transmural map clearly depicts BZ areas surrounding CZ regions (high probability areas in pink in Fig. 2f). During the intervention, potential ablation targets obtained from DE-MRI are analyzed by comparing them with tissue characterization (HT, CZ and BZ) performed on the corresponding EAM (BIP maps). For EAM the following standard thresholds are defined based on voltage criteria, as in (Andreu et al., 2011):

HT ≥ 1.5mV > BZ ≥ 0.5mV > CZ.

(3)

Using these standard thresholds, the average scar area (CZ + BZ) from EAM data in our patient population was of 59.74 ± 40.6 cm2 , in agreement with similar studies (Fernández-Armenta et al., 2013; Berruezo et al., 2015), which reported 57.6 ± 37.4 cm2 and 52 ± 39 cm2 , respectively. 3.2. Generation of synthetic ground-truth data The LV geometry of a randomly chosen patient from the available clinical dataset was taken as a reference surface to generate the ground-truth data. Two different scar configurations were simulated: (a) a scar located in the basal septum (top row in Fig. 3); and (b) a scar close to the apex (top row in Fig. 3). Both scar configurations included a CC and a small patch of BZ (non-CC) close to the CZ limits. Each scar was defined as a set of mesh points within

D. Soto-Iglesias et al. / Medical Image Analysis 32 (2016) 131–144

135

Fig. 3. Synthetic datasets. Top row: septal scar; bottom row: apical scar. Colors correspond to different tissue regions: core zone (CZ, red), border zone (BZ, green) and healthy tissue (HT, pink). n: number of points; σ : standard deviation of applied noise. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

a 1.5 cm radius circle (Euclidean distance in 3D space), representing an area of 7.78 cm2 (including a CZ area of around 6.42 cm2 and a BZ patch area 1.36 cm2 ) of the total LV surface area (19.20 cm2 ). The DE-MRI surfaces with synthetic scars were subsequently transformed four different times applying a random rigid transformation with a translation of up to 10 mm and rotation of up to ± 1 rad, as in Roujol et al. (2012), resulting in a set of rigidly transformed DE-MRI surfaces. From these surfaces, EAM data was simulated ensuring that the density of points and its local spatial variation (induced by the added noise) resemble real EAM characteristics. In order to model different densities of the EAM data, the synthetic DE-MRI surfaces were randomly downsampled using 300, 40 0, 50 0, 60 0 vertices, making sure that 50% of them were sampled at the scar and CCs (point density is usually higher in these zones). For simulating uncertainties in EAM point localization and imperfect catheter contact a Gaussian noise (standard deviation of σ = 0, 1, 2 and 3), as in Roujol et al. (2012), was also added to the non-scar vertex coordinates along the surface normal, inducing a maximum point displacement in the order of 9 mm. Due to the high density of points inside the scar a correlated Gaussian noise was added to the scar vertex coordinates along the surface normal. To impose noise correlation, 30 points within the scar were randomly selected and a random displacement was propagated to its 4 mm radius neighborhood as follows: d2

Displacement = Dmax · e(− 2·σ 2 )

(4)

where d is the distance of each point of the neighborhood to its center, σ = 0.5, and Dmax is the maximum displacement (along the surface normal) allowed. Maximum displacements of Dmax = ±1, ±2, ±3mm were tested. Each scar configuration was generated four times with different amounts of random noise and results were averaged in order to reduce the bias towards randomly created favorable or unfavorable displacements. Examples of DE-MRI surfaces with synthetic scars and simulated EAM ground-truth data can be seen in Fig. 3. Overall, a set of 128 synthetic datasets with different scar configurations (two), rigid transformations (four), EAM densities (four) and EAM point uncertainty (four) was generated to test the available registration and mapping strategies. All synthetic ground-truth datasets can be downloaded from the following website: http://physense.upf.edu. Besides studying the influence of EAM data acquisition parameters on the obtained integration results, the synthetically generated ground-truth data was used to perform a sensitivity analysis on the landmarks guiding the registration stage. The analysis was carried out by adding random noise of up to 1, 2 and 3 mm in the XYZ di-

rections to the landmark positions in the high-resolution DE-MRI mesh (10,0 0 0 vertices). 3.3. Accuracy measures The accuracy of the evaluated integration strategies was analyzed using: (i) global measures, providing a single figure-of-merit of the overall performance averaged across the different tissue labels; and (ii) local measures that allow a more specific regional analysis, especially in the scar area. The point-to-surface average distance between the registered EAM data points and the DE-MRI surfaces gives an initial estimation of the global accuracy of the different registration strategies. For the whole integration pipeline, including mapping algorithms, the Tanimoto Coefficient with Multiple Fractional labels (Crum et al., 2006), TCMF , an overlap measure that considers the presence of multiple heterogeneous and fuzzy labels, is well-suited to our problem since there are three different regions contributing to the overall accuracy (CZ, BZ and HT). It can be computed as follows:



T CMF



αl voxels,i MIN (Ali , Bli )  =  , α l abel s,l l voxels,i MAX (Ali , Bli ) l abel s,l

(5)

where α l is a label-specific weighting factor that affects how much each label contributes to the overlap accumulated over all labels. Ali and Bli are the amount of labels A and B at the voxel i. The relative area of each region over the total area of the LV surface mesh is estimated to obtain each α l as one minus the relative area for each region, thus giving enough weight to small regions (e.g. BZ). The accuracy of the different integration strategies can also be evaluated by studying local changes undergone around the scar. During the intervention double potential (DP) points are manually annotated after a careful visual analysis of the 1D electrograms in scar regions since they are usually located inside the scar and indicate the presence of a CC. In addition, entrances of CCs are also annotated in both the EAM and DE-MRI data. The following local accuracy measures were then compared: (a) the probability of each DP to belong to each region (CZ, BZ or HT); (b) the preservation of the relative area of each region after the integration; (c) the geodesic distance between the CC entrances; and (d) the CC entrance location error. 4. EAM and DE-MRI data integration with a quasi-conformal mapping (QCM) DE-MRI and LV EAM information from both modalities needs to be mapped onto the same reference space for comparison. It must

136

D. Soto-Iglesias et al. / Medical Image Analysis 32 (2016) 131–144

be pointed out that DE-MRI meshes are denser than EAM ones (number of vertices is 9230 ± 513.17 and 2204 ± 528.83 for DEMRI and EAM surface meshes, respectively). In addition, surfaces extracted from DE-MRI provide a more accurate anatomical representation of the LV compared to EAM data, which can be sparse and noisy. For these reasons, the DE-MRI LV surface is taken as reference space where to map the available information. Left ventricle surface meshes have a form of a cut ellipsoid and are homeomorphic to a 2D disk, i.e. there exists a continuous function between the LV surface and the disk with a continuous inverse (i.e. a bijective mapping) that preserves all topological properties. This homeomorphism can be calculated by requiring that every vertex coordinate of the triangulated surface mesh has a vanishing Laplacian, i.e. minimizing the difference between its position and the average (e.g. barycenter) of its neighborhood. By mapping both DE-MRI and EAM surfaces to the same disk we can then establish a homeomorphism between the two surfaces (De Craene et al., 2012). A disk representation of the LV is quite common in cardiology, often referred as to bull’s eye plot (Duchateau et al., 2013). Let ϕ : M → D ⊂ R2 be the sought bijective mapping between the surface M and a 2D unit disk D, being ∂ M and ∂ D their boundaries, respectively. In the particular case of LV surfaces, ∂ M coincides with the Mitral Annulus (MA). For simplicity, and without loss of generality, let us assume that D is in the XY plane. The following development is detailed for the x-coordinate but an equivalent procedure is also required for the y-coordinate to obtain for every vertex in the surface M its corresponding position in the disk D. Notice that the x-coordinates of ∂ D are given by a vectorcolumn x0 (concatenation of x coordinates of all boundary points) and that its size is defined by the mesh resolution of ∂ M for dimension consistency, i.e. it is equal to the number of points in ∂ M. Let LMࢨ∂ M be the Laplacian matrix of the surface mesh M without the rows corresponding to its boundary ∂ M. Let x∂ D and xDࢨ∂ D be the x-coordinates of the boundary and the interior points of the disk, respectively. Then, the solution of the following system of linear equations gives the sought bijective mapping ϕ :



LM\∂ M · x D\∂ D = 0 x∂ D = x0

(6)

The bijective mapping ϕ is defined by vertex correspondences provided by Eq. (6) and the connectivity retained from the mesh representation of M (Tutte, 1963), considering ϕ as a piecewise linear map based on the barycentric coordinates of triangles. This methodology provides a simple method for mapping surfaces with one hole to a unit disk, where the boundary of the surface is mapped to the circumference of the disk (x∂ D = x0 ), with the imposed constraint of preserving the relative distance between the points of the boundary. This relative distance is defined as the percentage of the length on every edge over the perimeter of the whole mesh boundary. For the Laplacian computation the approximation is based on cotangent weights (Pinkall and Polthier, 1993), which gives rise to a conformal map, i.e. a function which preserves angles locally. However, since after the mapping there is no guarantee that the cardiac apex will end up in the same location on the disk for different meshes (due to different triangulations), its location is enforced to be at the center of the disk using Thin-Plate splines (TPS). TPS provide smooth transformations while imposing some landmark-based constraints. In this work the disk boundary is fixed and the apical landmark is displaced to the centre of the disk, then TPS are used as an interpolating function to apply the resulting transformation to all vertices of the mesh. The use of TPS in the proposed methodology relaxes the conformal mapping onto a quasi-conformal mapping, as detailed in the Appendix, which is still convenient for our application.

The mappings of the DE-MRI and EAM surfaces to the disk (ϕ demri → D and ϕ eam → D , respectively) can be considered two dif−1 ferent triangulations of the same disk. At this point, ϕdemri · →D ϕeam→D can be used to map any scalar field from the EAM surface to the DE-MRI surface. The procedure is as follows. Let xeam → D and xdemri → D be the (x,y) coordinates of the vertices mapped onto the disk D for the EAM and DE-MRI surfaces, respectively. Then, the matrix T expressing every point xdemri → D in terms of barycentric coordinates of the triangulation of xeam → D is constructed based on the mapping ϕ . Let ceam be the scalar values associated to the vertices of EAM mesh (e.g. LAT or BIP values). The corresponding scalars on the DE-MRI mesh can then be computed as:

ceam→demri = T · ceam .

(7)

In the mapping obtained with the methodology described above there is an orientation ambiguity due to the symmetry of the LV. In order to produce an univocal mapping, a landmark is selected on the mesh boundary, ∂ M, for both surfaces. In this work the landmark is chosen in the middle of the basal septum (MS), which is then mapped to the point (-1,0) in the disk ∂ D (see Fig. 4). 5. State-of-the-art EAM and imaging integration strategies The developed evaluation framework has been applied to several state-of-the-art integration strategies to compare their performance with the one provided by the QCM-based technique. These different integration strategies have been tested on the available clinical data and the synthetically generated ground-truth dataset, giving results that are analyzed with the proposed local and global measures. Integration strategies usually follow a two-stage approach: (i) a registration between EAM and DE-MRI meshes; and (ii) a mapping of EAM information over the DE-MRI geometry. The following subsections explain these stages in more detail. 5.1. Registration techniques The registration process between EAM and DE-MRI data aims at spatially aligning them since LV geometries from these two modalities are in different systems of coordinates. In clinical practice, only rigid (including translations and rotations) registration methods, based on manually identified landmarks, are typically used. In this work, three different registration techniques with different complexity transformation models were studied: landmarkbased rigid registration; the Iterative Closest Point (ICP) (Besl and McKay, 1992); and a non-rigid registration technique based on currents (Vaillant and Glaunes, 2005). 5.1.1. Landmark-based registration It estimates a rigid transformation between the EAM and DEMRI endocardial surfaces based on a set of landmarks manually placed at corresponding points. In our case, seven landmarks are selected in both surfaces, six on the mitral annulus (midseptal, infero-septal, infero-lateral, midlateral, antero-lateral and anteroseptal) and an additional one at the apex, as can be seen in Fig. 4. 5.1.2. Iterative closest point (ICP) registration The landmark-based registration is quite sensitive to the choice of landmarks and their position in both LV surfaces. The Iterative Closest Point (ICP) (Besl and McKay, 1992) algorithm can be consequently applied to partially correct possible landmark misplacements. Therefore, ICP was applied after obtaining an initial transformation with the landmark-based registration.

D. Soto-Iglesias et al. / Medical Image Analysis 32 (2016) 131–144

137

Fig. 4. Manually selected landmarks in DE-MRI (left column) and EAM (right column) data guiding the registration and mapping processes in a clinical case: midseptal (MS), infero-septal (IS), infero-lateral (IL), midlateral (ML), antero-lateral (AL) and antero-septal (AS) parts of the mitral annulus; and one additional landmark at the apex. Top: 3D mesh visualization. Bottom: disk representation. The DE-MRI data colormap represents pixel intensities (PI, low and high contrast being red and blue, respectively) while the EAM one shows bipolar voltages (BIP, low and high voltages being red and pink, respectively). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

5.1.3. Non-rigid registration based on currents To complement the previous two registration methods, a non-rigid registration technique based on currents (Vaillant and Glaunes, 2005) was also studied. This technique aims at minimizing the distance between two surfaces in the Large Deformation Diffeomorphic Metric Mapping framework (LDDMM) (Beg et al., 2005). This technique provides diffeomorphic transformations preserving the topology of the initial surfaces (in our application, the LV surfaces) and does not require point correspondences between the EAM and DE-MRI LV surfaces. Either landmark- or ICPbased rigid registration techniques can be used to provide an initial transformation to the LDDMM algorithm. 5.2. Mapping techniques We evaluated the influence of the mapping technique on final integration results comparing the most commonly used mapping technique (closest point, CP) with the new QCM-based one proposed in Section 4. Given the two registered surfaces represented by triangular meshes, the BIP value of every vertex of the DE-MRI mesh needs to be assigned. Every vertex of the DE-MRI mesh is projected to the closest face over the EAM mesh. Then, the BIP values associated to the vertices of this EAM face are linearly interpolated and associated to each original DE-MRI vertex. 6. Results In this section, results obtained from applying the evaluation framework on the different integration strategies, including the synthetic ground-truth (128 datasets) and clinical (17 cases) data of VT patients are presented.

Table 1 Point-to-surface errors provided by the different registration techniques applied to synthetic datasets. Bold values indicate the registration method providing best results. Point-to-surface error

Landmarks ICP Currents

Apical scar

Septal scar

2.02 ± 1.87 mm 1.24 ± 0.90 mm 0.66 ± 0.26 mm

1.38 ± 0.73 mm 1.04 ± 0.66 mm 0.65 ± 0.24 mm

6.1. Synthetic ground-truth data 6.1.1. Accuracy of registration techniques Registration accuracy from the three tested techniques was evaluated on the whole set of synthetic data computing the pointto-surface average distance between the EAM (points) and the DEMRI (surfaces) datasets, as reported in Table 1. Average errors are within an acceptable range for all tested registration techniques (mostly below 2 mm) considering EAM measurement uncertainties on spatial localization and imperfect catheter contact (they can easily reach 5 mm displacements). The best results were obtained with the current-based (non-rigid) registration technique that provided sub-millimeter accuracy, while the landmark-based registration produced larger errors in the apical scar configuration probably due to the lack of landmarks in this area. A Wilcoxon signedrank test was used on the point-to-surface error results to estimate the significance of differences between the registration methods, with a p-value <0.05 to determine statistically significant differences. The three registration methods achieved statistically different results.

138

D. Soto-Iglesias et al. / Medical Image Analysis 32 (2016) 131–144

Fig. 5. Multiple fractional generalized Tanimoto (mean ± standard deviation)) provided by applying the different integration methodologies on the synthetic datasets. The quasi-conformal (QCM) and closest point (CP) mapping techniques are combined with landmarks-based (Land), Iterative Closest Point (ICP) and currents registration algorithms to generate different integration methodologies. Results have been separated for different levels of noise (σ = 0, σ = 1, σ = 2 and σ = 3 in different color columns) and for different scar configurations (A, septal scar; B, apical scar). Table 2 Local evaluation measures obtained with different integration methodologies applied to synthetic ground-truth datasets with different noise levels: landmarks (Land), Iterative Closest Point (ICP) and currents registration techniques with closest point (CP) or quasi-conformal (QCM) mapping stages. CZ: core zone; BZ: border zone; HT: healthy tissue; CC: conduction channel. Bold values indicate the best result among the different integration alternatives. Noise

Area preservation (%) CZ

Land-CP

ICP-CP

Currents-CP

Land-QCM

Currents-QCM

0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3

16.66 19.37 17.12 21.17 17.96 19.05 20.01 23.43 60.33 62.20 61.83 65.74 −5.96 −5.50 −3.57 −1.16 37.34 39.35 38.96 42.78

BZ ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

9.43 10.19 10.63 16.11 8.63 8.68 9.84 12.77 22.33 22.08 23.92 22.78 7.82 8.44 10.86 12.67 18.74 17.55 18.13 18.28

62.63 59.24 59.76 46.69 65.11 59.23 55.70 49.77 123.46 121.09 113.90 114.41 18.34 19.39 17.48 17.58 75.20 78.30 72.48 75.72

HT ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

30.90 32.13 38.19 36.93 29.25 32.22 33.96 34.87 65.63 64.70 69.59 67.29 12.13 11.63 11.75 12.60 49.15 44.73 45.62 43.99

6.1.2. Accuracy of the whole integration pipeline The Multiple Fractional labels version of the Tanimoto Coefficient (TCMF ) was used here to assess the global accuracy of the different evaluated integration pipelines, which is summarized in Fig. 5. In this figure, the TCMF is represented as the mean ± standard deviation averaged over all simulated EAM data (different rigid transformations and point densities) for the two scar configurations (septal and apical scars on the left and right of Fig. 5, respectively) and the four different amounts of noise added to the synthetic datasets (different columns in Fig. 5). The different combinations of registration and mapping techniques tested on the synthetic datasets were: landmark-based registration with both closest point and quasi-conformal mappings (Land-CP and Land-QCM, respectively); ICP registration with closest point mapping (ICP-CP)2 ; and currents-based non-rigid registration with both closest point and quasi-conformal mappings (Currents-CP and Currents-QCM, respectively).

2 ICP-QCM and Land-QCM results are equivalent since the QCM technique makes integration only depending on the LV geometry and the landmark selection not on the spatial position of the mesh (up to rigid transformations).

−0.86 −0.91 −0.86 −0.88 −0.91 −0.91 −0.91 −0.96 −2.44 −2.48 −2.42 −2.53 0.04 0.02 −0.02 −0.08 −1.51 −1.58 −1.53 −1.66

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.29 0.28 0.39 0.49 0.24 0.24 0.26 0.31 0.95 0.93 1.01 0.89 0.17 0.19 0.27 0.35 0.79 0.72 0.74 0.69

CC entrance location

CC length

error (mm)

difference (mm)

4.06 4.23 4.40 4.66 3.80 3.78 4.11 4.30 7.20 7.18 8.04 7.23 1.83 2.45 2.59 3.31 4.67 5.42 5.28 5.62

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

2.12 2.43 2.23 2.53 2.17 1.93 2.12 2.34 4.75 4.21 5.75 3.84 1.50 1.68 1.72 1.98 4.29 4.31 3.63 3.46

4.53 4.96 4.52 4.54 4.72 4.66 4.82 4.88 9.13 9.33 9.37 8.96 1.06 1.52 1.29 1.63 5.82 6.49 6.07 6.00

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

2.07 2.00 2.33 2.42 2.06 1.93 1.92 1.85 4.87 4.32 4.24 3.80 1.35 1.52 1.37 1.81 4.43 4.01 3.66 3.60

Fig. 5 shows that the Land-QCM consistently provided the best global accuracy of the integration accuracy in every experiment with synthetic datasets. A Wilcoxon signed-rank test was used on the TCMF results to estimate the significance of differences between the proposed integration pipeline (Land-QCM) and the other alternatives, with a p-value < 0.05 to determine statistically significant differences. Results obtained with the Land-QCM pipeline were statistically different from all other methods. It can also be observed in Fig. 5 that the QCM mapping technique outperforms the CP one independently of the registration method, demonstrating an improvement in accuracy with both landmarks- and currents-based registration techniques. Another interesting observation from these results is that the use of complex non-rigid registration techniques does not improve the integration accuracy compared to the simple rigid landmark-based registration method when concatenating with the QCM mapping approach. Finally, there are not substantial differences in overall performance of the integration techniques between the two different scar configurations. Local measures of accuracy at the scar level are summarized in Table 2. Similar to the global accuracy analysis described above

D. Soto-Iglesias et al. / Medical Image Analysis 32 (2016) 131–144

139

Fig. 6. Integration of EAM and DE-MRI data after applying the different evaluated strategies on synthetic ground-truth data (two left columns) with two scar configurations (septal and apical scar, top and bottom rows, respectively), 300 EAM simulated points and noise of σ =3. Colors correspond to different tissue regions: core zone(CZ) in red, border zone (BZ) in green and healthy tissue (HT) in pink. Landmarks (Land), Iterative Closest Point (ICP) and currents registration techniques combined with closest point (CP) or quasi-conformal mapping (QCM) mapping methods provided the evaluated integration methodologies. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the proposed integration methodology, Land-QCM, again outperforms any other combination of registration and mapping techniques for all local accuracy measures in the different myocardial tissues and for CC characterization. This is also achieved independently of the amount of noise added to the data. Results obtained with the Land-QCM pipeline were statistically different (p-value < 0.05) from all other methods for area difference, CClength difference and CC-entrance location error criteria. The improvement is particularly striking when comparing with the current standard way of integrating DE-MRI and EAM data in clinical routine (Land-CP or ICP-CP strategies). For instance, area differences in the BZ are 65.69% ± 34.27 vs 56.62% ± 31.37 vs 17.48% ± 13.25 for Land-CP, ICP-CP and Land-QCM, respectively. Similarly to observations made on global accuracy measures the non-rigid registration approach provided the worst results, partially corrected when switching from the CP to the QCM mapping method. Fig. 6 provides visual support to the quantitative results provided by global and local measures of accuracy on the performance of the evaluated integration strategies. In this figure, an illustrative example of the synthetically generated ground-truth data, including the DE-MRI LV surface with synthetic scars and a simulated EAM with 300 points and noise σ = 3 is shown for the two scar configurations (septal scar, top row; apical scar, bottom row; Fig. 6). Integration results provided by the different evaluated approaches are superimposed on the common reference space (DEMRI surface) for comparison purposes. Several observations can be made from a visual inspection of these results. For instance, scar volume is incorrectly larger when applying the CP rather than with the QCM mapping technique. Additionally, one can see in the figure that CC main characteristics (connectivity, shape, location, extent) are better preserved with the Land-QCM integration technique. The non-rigid registration method deforms the EAM mesh basically using mesh point location information. Therefore, most spatially different corresponding EAM and DE-MRI points will guide the estimation of the transformation model, independently on the tissue label of each vertex, which makes scar tissue suffer unrealistic large deformations if differences between EAM and DE-MRI surfaces are large around the scar (note the abnormally larger CCs in currents-based approaches in the last two columns of Fig. 6). 6.1.3. Influence of landmark selection Results from the study of the influence of landmark selection errors on the integration pipeline performance are presented in

Table 3 Generalized Tanimoto coefficient obtained after applying the evaluated integration methodologies on synthetic data with different landmark selection errors: landmarks (Land), Iterative Closest Point (ICP) and currents-based registration techniques combined with closest point (CP) or quasi-conformal mapping (QCM) mapping methods. Apical scar

Land-CP

ICP-CP

Currents-CP

Land-QCM

Currents-QCM

Septal scar

Error (mm)

Septal Landmark

Apical Landmark

Septal Landmark

Apical Landmark

1 2 3 1 2 3 1 2 3 1 2 3 1 2 3

0.97 0.90 0.85 1 0.88 0.85 1 0.90 0.86 0.98 0.94 0.90 0.99 0.97 0.97

0.90 0.88 0.84 1 1 1 1 1 1 0.89 0.87 0.82 0.99 0.99 0.99

0.93 0.83 0.76 1 0.82 0.77 1 0.83 0.78 0.89 0.80 0.71 1 0.94 0.90

0.97 0.94 0.91 1 1 1 1 1 1 1 1 0.99 1 1 1

Table 3. The TCMF coefficient for measuring global accuracy was used here. It can be observed in the table that the integration performance gets worst with larger landmark errors, as expected. Moreover, these findings suggest that the correction of these landmark errors are mainly associated to the accuracy of the registration stage rather than the mapping one since the non-rigid registration strategies perform substantially better than the rigid transformation ones. Intuitively, landmark-based registration methods are the most sensitive to landmark positioning errors, as confirmed in Table 3. It can also be observed that rigid registration with ICP can recover landmark errors up to a level of noise (in general, to 2 standard deviations), but with larger landmark errors, non-rigid registration is required. Comparing both CP and QCM mapping methods, Land-QCM is not as robust as Land-CP to landmark errors in some particular cases (e.g. apical scar and apical landmark errors). This is due to the Land-QCM pipeline, which uses the apical landmark for the recentering of the apex, making it more dependent to its selection. In addition, septal landmark errors are generally more critical than apical ones.

140

D. Soto-Iglesias et al. / Medical Image Analysis 32 (2016) 131–144 Table 4 Generalized Tanimoto coefficient, regional area differences (mean% ± STD%) and probability of double potential (DP) EAM points belonging to scar (BZ and CZ) or HT regions on the clinical cases for the different integration methodologies: landmarks (Land), Iterative Closest Point (ICP) and currents registration techniques combined with closest point (CP) or quasi-conformal mapping (QCM) mapping methods. CZ: core zone; BZ: border zone; HT: healthy tissue. TCMF

Area differences (%) CZ

Land-CP ICP-CP Currents-CP Land-QCM Currents-QCM

0.22 0.23 0.24 0.27 0.25

± ± ± ± ±

0.09 0.09 0.09 0.09 0.09

28.74 34.61 32.78 27.12 29.34

Probability of DP points in each tissue label (%)

BZ ± ± ± ± ±

27.22 34.60 26.57 19.72 21.45

22.12 16.32 16.54 16.24 18.66

HT ± ± ± ± ±

24.35 15.60 17.31 19.76 19.08

8.21 6.82 7.14 6.15 7.94

6.2. Clinical data 6.2.1. Accuracy of registration techniques The point-to-surface average distance between the EAM (points) and DE-MRI (surfaces) data for clinical cases registered with the techniques described above were the following: (i) Landmarks, 5.45 ± 1.21 mm; (ii) ICP, 3.60 ± 0.86 mm; (iii) Currents, 0.74 ± 0.20 mm. As with synthetic ground-truth data, the nonrigid registration technique based on currents achieved the best fitting results, even though the overall performance for all registration techniques is worst on clinical data than in synthetic datasets. It must also be pointed out that registration errors obtained with rigid transformation models are of the same order as published in (Tao et al., 2012; Roujol et al., 2012) (average errors of around 5 mm). 6.2.2. Accuracy of the integration pipeline Some of the global and local evaluation measures applied to synthetic ground-truth data can be used on clinical data but they need to be carefully interpreted and considered as a semiquantitative estimation of the integration performance. The global accuracy of the integration approaches applied on clinical data was estimated with the TCMF between original EAM and DE-MRI values and by estimating regional area preservation after tissue classification based on the standard clinical thresholds in both modalities. For each region we computed its area relative to the whole LV surface area, before and after applying the different integration pipelines in the EAM data. Changes in regional relative areas indicates how well tissues are preserved after integration with DE-MRI data. The obtained results are summarized in Table 4. In general the proposed integration methodology (Land-QCM) performed better than the other techniques but differences are not as important as when processing the synthetic datasets. Nevertheless, we can observe a clear improvement in the TCMF coefficient of the QCM over the CP mapping techniques when applied to clinical data. It should also be noted that the TCMF values obtained here are in the same range to those obtained when comparing scar segmentation techniques with manually obtained ground-truth in the DE-MRI challenge recently organized at STACOM’123 , giving an order on the maximum expected accuracy with the current data resolution when analyzing these patchy scar structures. Additionally, relative area differences for each region after threshold-based tissue classification show a general superior performance of the Land-QCM strategy, achieving better HT and CZ preservation. However, the BZ preservation is in the same range than ICP-CP, being the strategies that better preserve the BZ relative area in our experiments (16.24 ± 19.76 and 16.32 ± 15.60 percentage of area differences for Land-QCM and ICP-CP techniques, respectively). The BZ preservation is quite important in our application since it often contains small structures within the scar such as CCS that will be used to guide the ablation intervention. The application of non-

3

http://www.physense.org/stacom2012.

scar ± ± ± ± ±

3.81 3.23 4.52 5.31 4.41

67.62 72.95 73.47 74.33 73.40

CZ ± ± ± ± ±

19.30 16.68 15.73 18.47 19.91

28.00 32.21 32.29 28.97 30.85

BZ ± ± ± ± ±

23.18 22.76 20.20 19.99 22.16

39.52 40.65 41.26 45.35 42.54

HT ± ± ± ± ±

13.40 13.09 14.30 17.61 16.58

32.48 27.06 26.53 25.67 26.60

± ± ± ± ±

19.30 16.68 15.73 18.47 19.91

rigid registration methods (e.g. currents-based) tends to incorrectly change the relative proportion of HT with respect to scar areas independently of the mapping method applied, as shown by the increment of HT preservation errors shown in Table 4: from 6.82 ± 3.23 (ICP-CP) to 7.14 ± 4.52 (currents-CP); and from 6.15 ± 5.31 (Land-QCM) to 7.94 ± 4.41 (currents-QCM). An additional local measurement was also computed to estimate the performance of the integration techniques in clinical cases with CCs identified in EAM data. The location of EAM DP points (blue dots in Fig. 7) in the corresponding DE-MRI surface after applying the integration procedure was analyzed depending on the underlying tissue transmural map (HT, BZ, CZ; but also including a scar label combining CZ and BZ regions). The average number of DP points per case was 54.01 ± 30.32 and the obtained results are shown in Table 4. It can be observed that the QCM mapping technique performs better than the CP one in all regions, locating DP points mainly in the scar region defined by the DE-MRI data (74.33% of DP points). However, identifying the best result within the scar is not straightforward since some DPs should correctly be located in CZ rather than BZ regions, due to the current DE-MRI resolution limitations (i.e. some small CCs cannot be seen). In any case the QCM-based integration techniques locate more DP points in BZ areas than the CP ones, being the Land-QCM the one providing the highest probability. Figs. 7 and 8 give two illustrative examples on the integration results with clinical data. It can be appreciated the difficulties of evaluating the integration accuracy for some cases where DE-MRI and EAM data, especially scar information, do not look very similar. There could be several reasons for these differences such as the inappropriateness of standard threshold values or different fields of view of the integrated data. Fig. 7 depicts data from a clinical case with a scar presenting a small CC, which is not visible in DEMRI, while Fig. 8 shows a scar where the CC is large enough to be recognized in DE-MRI. Results of the integration strategies shown in these figures demonstrate that QCM better preserves EAM scar characteristics (area, CCs) in both clinical cases, whereas the CP mapping can abnormally modify them, depending on the LV surface differences between the EAM and the DE-MRI data. A visual inspection of the manually annotated DP points (blue and black dots in the figures) after applying the integration pipelines points out that the CP mapping technique consistently misplaces EAM points outside the DE-MRI CC, especially when combined with rigid registration techniques. For instance, we can observe in Fig. 7 two mitral annulus points (white dots, numbers 1 and 2, top part of the LV surface) in the EAM data, which are incorrectly located in the integrated DE-MRI surface when using CP combined with rigid registration techniques while the proposed integration method (Land-QCM) successfully places them in the annulus. The Land-QCM also succeeds at locating most of the EAM DP points within the CC defined with DE-MRI. The manually identified entrances of the CC from the EAM data (black points 3 and 4) are mostly located closer to the DE-MRI ones using the LandQCM strategy than with the other integration alternatives. It must

D. Soto-Iglesias et al. / Medical Image Analysis 32 (2016) 131–144

141

Fig. 7. Integration of EAM and DE-MRI data from a clinical case with a simple scar having a small conduction channel (CC, white line) only identified in EAM. Top row: (a) EAM data visualized in the CARTO viewer; (b) EAM data with tissue characterization obtained with standard thresholds (core zone, CZ in red; border zone, BZ, in green; healthy tissue, HT, in pink); (c–g) EAM tissue labels mapped on the DE-MRI geometry after applying the different integration methodologies (quasi-conformal, QCM, and closest point, CP, mapping techniques are combined with landmarks-based, Land, Iterative Closest Point, ICP, and currents registration algorithms). Middle row: (a) 10% layer of the DE-MRI data with tissue characterization obtained with standard thresholds; (b) CZ transmural map, from 0 (pink) to 1 (red); (c–g) Annotated EAM points mapped over DE-MRI geometry after applying the different integration methodologies. Double potential points (green), ablation points (red), CC entrances (black, e.g. points 3 and 4), mitral annulus points (white, e.g. points 1 and 2) are shown. Bottom row: zoom of middle row around the conduction channels.

be pointed out that this CC is relatively small and subsequently can only be detected after careful analysis of the 1D electrograms of EAM points around this area, but it cannot be identified in DEMRI data. This type of situations could explain the fact that a high number of DP points end up in the CZ region (rather than in BZ) after applying Land-QCM (see Table 4). Fig. 8 shows an example where EAM and DE-MRI meshes have different fields of view since basal planes in both modalities are at different levels, being the EAM mesh size smaller than the DEMRI one. It can be seen in the figure how the CP mapping cannot correct these differences, especially when combined with rigid registration techniques, and tends to accumulate EAM points in areas where the point to surface error is larger, thus wrongly projecting several EAM points to the same closest DE-MRI location (see cloud of points in the mid basal region, 1 in Fig. 8). This accumulation of points is partially reduced after applying non-rigid registration due to the reduction of the point-to-surface error between the EAM and DE-MRI LV geometries. Analyzing the results around the CC, we can observe that the Land-QCM method mostly places DP points where the CC is defined in DE-MRI data, whereas the remaining integration pipelines fail, incorrectly moving them to the apex (ICP-CP, Currents-CP, Currents-QCM) or to a more basal position (Land-CP). When focusing on the CC entrances (black points 2, 3 and 4 in Fig. 8), which are potential targets for ablation, it can be seen how Land-QCM locates them in the extremities around the CC, where the remaining methods generally place them farther from the CC. 7. Discussion and conclusions There is a high complementarity in myocardial tissue information acquired from EAM and DE-MRI modalities, as it is clearly visible in Figs. 7 and 8, since they represent different biophysical phenomena that are studied at different times and conditions. Other

differences between the two modalities also appear due to the use of different coordinate and resolution systems, which can be partially compensated using computational techniques to spatially integrate both types of data. Nevertheless, performing a rigorous validation of the accuracy of these integration techniques is not a trivial task due to the biophysical differences between the modalities that do not need to be compensated and the lack of ground-truth data to assess in a quantitative way the performance of the whole computational process. In this paper, we present an evaluation framework that allows a rigorous comparison of different EAM and DE-MRI integration strategies, also being able to identify which stages are more critical for the final results. Clinical and synthetically generated ground-truth data in combination with a set of global and local accuracy measures were used to assess state-of-the-art integration techniques and the QCM-based one proposed in this paper. The developed synthetic dataset of 128 different simulated EAMs, together with the corresponding MRI data, permitted a fully quantitative evaluation of the ability of the integration techniques to cope with spatial differences between EAM and DE-MRI coordinate systems, neglecting other biophysically-related differences present in clinical data. This ground-truth dataset can be further developed to consider more realistic scar configurations and then be used for benchmarking purposes with more advanced integration techniques. Global accuracy measures used in the evaluation framework were quite useful to give an overview of the performance of each integration technique but in our application a regional analysis of the integration accuracy is required. For this reason we developed a set of local measures that successfully relate EAM critical information (e.g. double potential points) with scar information extracted from DE-MRI (scar transmurality). They allowed a local analysis of the integration accuracy in the regions of interest, i.e. basically in core and border zone regions to assess the misplacement of critical targets for ablation therapy.

142

D. Soto-Iglesias et al. / Medical Image Analysis 32 (2016) 131–144

Fig. 8. Integration of EAM and DE-MRI data from a clinical case with a complex scar having a conduction channel (CC, white line) visible in DE-MRI. Top row: (a) EAM data visualized in the CARTO viewer; (b) EAM data with tissue characterization obtained with standard thresholds (core zone, CZ in red; border zone, BZ, in green; healthy tissue, HT, in pink); (c–g) EAM tissue labels mapped on the DE-MRI geometry after applying the different integration methodologies (quasi-conformal, QCM, and closest point, CP, mapping techniques are combined with landmarks-based, Land, Iterative Closest Point, ICP, and currents registration algorithms). Middle row: (a) 10% layer of the DE-MRI data with tissue characterization obtained with standard thresholds; (b) CZ transmural map, from 0 (pink) to 1 (red); (c–g) Annotated EAM points mapped over DE-MRI geometry after applying the different integration methodologies. Double potential points (green, e.g. cloud of points 1), ablation points (red), CC entrances (black, e.g. points 2, 3 and 4) are shown. Bottom row: zoom of middle row around the conduction channels. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Until now, integration approaches have just focused on improving the registration stage, using the simplest closest point approach for the subsequent mapping phase. Using the developed evaluation framework we demonstrate in this work that the mapping stage is as relevant as the registration one for a successful integration. Several integration strategies were compared using different registration methods and mapping techniques. As can be seen in Table 2 and Fig. 5 the proposed integration pipeline (Land-QCM) provides the best results in all local and global measurements estimated from synthetic data (better area preservation, smaller error in entrance localization and CC length), reducing errors by a factor of 50-100% with respect to state-of-the-art techniques currently used in clinical routine. QCM preserves the topology of the mapped structures while the CP technique does not assure this property. This preservation is critical when mapping scars with conduction channels since the CP may make some CCs disappear. Furthermore, QCM is better than the typically used CP approach independently of the chosen registration technique. In addition QCM does not require complex registration algorithms since it is independent of the differences in coordinate systems of the data. Nevertheless, we would like to point out that the currents-based registration technique has been applied with default parameters and even though it already provides the best point-to-surface registration performance, it is likely that better results would have been obtained after some parameter tuning. The global overlap measures are still reasonable in synthetic data they are extremely difficult to interpret in clinical data. In the clinical dataset shown in Table 4, the Land-QCM still outperforms the alternative methodologies in terms of TCMF . It can also be seen that area differences highly depend on the evaluated region. Even if Land-QCM achieves the best results for all the regions, the BZ results are into the same range than CP-ICP and CPcurrents. This might be due to the EAM BZ topology, which is very different from DE-MRI one in some cases, as can be seen in Fig. 7. These differences in BZ can be explained by several factors such as the far field effect presented in the EAMs or the interpolation of

the EAM system if a small number of measurements are acquired. The analysis of the location of the DP points (manually identified in EAMs) in the DE-MRI surface after applying the different integration pipelines reveals that the proposed methodology correctly places them around the CC defined in the DE-MRI data. The other integration methodologies tend to move EAM points to wrong locations (e.g. outside the mitral annulus or the CC, as can be seen in Fig. 7), especially due to large point-to-surface errors when using rigid registration techniques. This issue is partially corrected with more complex non-rigid registration methods such as currents. The use of TPS to correct for apical misalignment converts QCM onto a quasi-conformal mapping since TPS will change small circles onto small ellipses but with bounded eccentricity (see Appendix). Other quasi-conformal mapping strategies have been recently published such T-map techniques proposed by Lui et al. (2014), which could also be used in our context. However we decided to use the TPS since the surface deformation induced by the TPS has a more global support than T-map ones, allowing for a better visualization of the data. T-map techniques tend to constrain mesh deformation to the smallest neighborhood of the displaced point, then compressing a lot of information in small areas, which is not convenient in our application. On the other hand, the main limitation of our methodology is its sensitivity to the septal landmark selection. As can be seen in Table 3 results provided by Land-QCM are influenced by landmark errors in some regions, such as in apical scars with apical landmark errors. However, with an appropriate scar characterization, these errors could be partially corrected by applying QCM after currents-based registration using scar information. The extension of the method to other cardiac chambers, such as RV and atria will also help on improving the robustness of the proposed methodology with respect to landmark errors. The proposed QCM technique establishes a common reference frame allowing the comparison of inter- and intra-patient data from different sources at different time-points. For instance, simulation results from computational models of cardiac electrophysiology, which can provide very useful clinical insight for op-

D. Soto-Iglesias et al. / Medical Image Analysis 32 (2016) 131–144

timizing RFA interventions (Relan et al., 2011; Smith et al., 2011; Ashikaga et al., 2013; Trayanova and Chang, 2016), could easily be compared with EAM and DE-MRI data in the common reference space. This would be beneficial both for making electrophysiological simulations more realistic but also to use them to compensate for possible data acquisition errors and then provide to clinicians better quality and more complete information for their clinical decisions. The accuracy of any integration technique is dependent on the quality of EAM and DE-MRI acquisitions, in particular in datasets with large differences in scar multi-modal information or fields of view, low spatial resolution or insufficient density of measurements. Nevertheless, there are continuous technological advances in EAM and MRI systems that are constantly reducing these issues. For instance, the use of multiple electrode catheters (e.g. Pentaray NAV catheter, BioSense Webster) together with better localization systems are allowing the acquisition of very dense EAMs (more than 10 0 0 points) in very short time (few minutes). Advanced signal processing techniques are also currently being developed to automatically process EGMs (Alcaine et al., 2014), which will help to detect outlier points and better characterize CZ and BZ regions. From the MRI side, the quality of 3T-scanners and the application of advanced MR image reconstruction algorithms, such as the one recently proposed by Ukwatta et al. (2015), significantly improves the 3D characterization of the scar area and the analysis of conduction channels since high-resolution isotropic images can be reconstructed (Andreu et al., 2015). These system advances are easing the task of finding spatial correspondences between both modalities, thus forcing state-of-the-art integration techniques to provide better accuracy. Tissue classification and scar segmentation will also be more robust in EAM and DE-MRI, which is not currently the case. In our experiments tissue characterization provided by standard clinical thresholds was not always optimal, leading us to perform a threshold scanning within a range of ± 10% around the standard clinical thresholds to visually maximize the overlap between the scar defined in EAM and DE-MRI, as suggested in Andreu et al. (2011). The lack of reliability of scar segmentations was the main reason for not using them to guide the integration techniques evaluated in this work. Roujol et al. (2012) proposed to constrain registration techniques guided by LV geometries with scar information but errors in tissue characterization are then propagated to the integration process. Nevertheless, with more robust scar segmentations this approach would be beneficial for the integration accuracy to avoid incorrect scar distortions when LV geometries from EAM and DE-MRI are very different (e.g. distinct definition of basal planes). From the other hand, a good integration strategy would also help to detect segmentation errors in both modalities, allowing a multi-modal approach for tissue characterization. In most clinical centers imaging data is mainly used in a qualitative way nowadays, for pre-operative planning purposes, searching for regions of interest with potential ablation candidates. Current limitations of MRI and EAM acquisition systems, together with the use of relatively simple image processing techniques, hamper a more extensive use of image-guided RFA interventions in VT cases, which require: (i) good quality MR images with high and (almost) isotropic resolution; (ii) strong experience with imaging technologies by electrophysiologists; (iii) advanced signal and image processing tools including tissue characterization and integration techniques, mainly in collaboration with engineering teams. Unfortunately, these requirements are only fulfilled in specialized and very active research units of arrhythmias. Nevertheless, in an appropriate environment an accurate integration technique can be applied to provide the clinician with MRI data fused with EAM information during the intervention itself. This integrated data is extremely useful intra-operatively to properly understand EAMs, which are

143

generally difficult to interpret without imaging information, and then having an impact on the decision on where to ablate. Ongoing randomized clinical studies at Hospital Clínic de Barcelona are giving promising preliminary results showing the clear benefit of integrating MRI and EAM data, significantly reducing the percentage of events after VT ablation. Technological advances in MRI and EAM systems will definitively make image-guided VT ablations possible in more clinical centers in the near future. The developed evaluation framework could also contribute in this issue by providing a pipeline for assessing the different steps of the data processing pipeline involved in image-guided VT ablations. The proposed QCM-based integration technique could also help on incorporating into the process in-silico indices derived from computational models of electrophysiology, among other types of complementary data. Acknowledgments This study was partially funded by Spanish Ministry of Science and Innovation (TIN2011-28067) and from the Seventh Framework Programme (FP7/2007-2013) for research, technological and demonstration under grant agreement VP2HF (no. 611823). Appendix A. Quasi-conformality of Thin-Plate Splines Thin-Plate Splines (TPS) provide a spatial mapping f : R2 → R2 that is defined as Bookstein (1989):

f (x, y ) = u(x, y ) + iv(x, y ) u(x, y ) = ax0 + ax1 x + ax2 y +



wxi ϕ (||Pi − (x, y )|| )

i

v(x, y ) = ay0 + ay1 x + ay2 y +



wyi ϕ (||Pi − (x, y )|| ),

(A.1)

i

where Pi are the coordinates of the control points, ai and wi are the weights (estimated during construction of the TPS from a set of point correspondences), ϕ (r ) = r 2 ln(r ) is the TPS radial basis kernel and ||·|| is the l2 norm. A mapping f is quasi-conformal if it satisfies the Laplace– Beltrami equation (Bers, 1977):

∂z¯ f = μ(z )∂z f,

(A.2)

with μ(z) being a Lebesgue measurable complex-valued function and supz |μ(z)| < 1, and z = x + yi. The derivatives can be defined through the Euler formulas:

∂z¯ f = 0.5([∂x u − ∂y v] + i[∂y u + ∂x v])

(A.3)

∂z f = 0.5([∂x u + ∂y v] + i[∂x v − ∂y u]).

(A.4)

By simple substitution it can be shown that

|μ ( z )|2 =

|∂z¯ f |2 ||J||2F − 2det (J ) = , |∂z f |2 ||J||2F + 2det (J )

(A.5)

where J is the Jacobian of the mapping f:

 ∂x u J= ∂x v

 ∂y u ∂y v

(A.6)

and ||J||F is the Frobenius norm. It can be seen that for |μ(z)| < 1 to hold for any z the mapping must be orientation preserving, i.e. det(J(z)) > 0 for every z. This implies that point correspondences defining the TPS must not have any change of orientation (flipping about an axis) nor surface folding (which would lead to loss of bijectivity). In the scenario treated in this paper where the boundary is convex and fixed and only one point is displaced within the boundary, the TPS will not produce any foldings. A more detailed treatment of the bijectivity ˚ of TPS can be found in Eriksson and Aström (2012).

144

D. Soto-Iglesias et al. / Medical Image Analysis 32 (2016) 131–144

References Alcaine, A., Soto-Iglesias, D., Calvo, M., Guiu, E., Andreu, D., Fernández-Armenta, J., Berruezo, A., Laguna, P., Camara, O., Martínez, J., 2014. A wavelet-based electrogram onset delineator for automatic ventricular activation mapping. IEEE Trans. Biomed. Eng. 61, 2830–2839. Andreu, D., Berruezo, A., Ortiz-Pérez, J., Silva, E., Mont, L., Borràs, R., de Caralt, T., Perea, R., Fernández-Armenta, J., Zeljko, H., Brugada, J., 2011. Integration of 3D electroanatomic maps and magnetic resonance scar characterization into the navigation system to guide ventricular tachycardia ablation. Circ.: Arrhythmia Electrophysiol. 4 (5), 674–683. Andreu, D., Ortiz-Pérez, J., Fernández-Armenta, J., Guiu, E., Acosta, J., Prat– González, S., Caralt, T.D., Perea, R., Garrido, C., Mont, L., Brugada, J., Berruezo, A., 2015. 3d delayed-enhanced magnetic resonance sequences improve conducting channel delineation prior to ventricular tachycardia ablation. Europace 17, 938–945. Ashikaga, H., Arevalo, H., Vadakkumpadan, F., Blake, R., III, Bayer, J., Nazarian, S., Zviman, M.M., Tandri, H., Berger, R., Calkins, H., Herzka, D., Trayanova, N., Halperin, H., 2013. Feasibility of image-based simulation to estimate ablation target in human ventricular arrhythmia. Heart Rhythm 10, 1109–1116. Auzias, G., Lefèvre, J., Troter, A.L., Fischer, C., perrot, M., Règis, J., Coulon., O., 2013. Model-driven harmonic parametrization of the cortical surface: HIP-HOP. IEEE Trans. Med. Imaging 32, 873–887. Beg, M.F., Miller, M.I., Trouvé, A., Younes, L., 2005. Computing large deformation metric mappings via geodesic flows of diffeomorphisms. Int. J. Comput. Vis. 61 (2), 139–157. Berruezo, A., Fernández-Armenta, J., Andreu, D., Penela, D., Herczku, C., Evertz, R., Cipolletta, L., Acosta, J., Borras, R., Arbelo, E., Tolosana, J., Brugada, J., Mont, L., 2015. Scar dechanneling: new method for scar-related left ventricular tachycardia substrate ablation. Circ.: Arrhythmia Electrophysiol. 8, 326–336. Berruezo, A., Fernández-Armenta, J., Mont, L., Zeljko, H., Andreu, D., Herczku, C., Boussy, T., Tolosana, J., Arbelo, E., Brugada, J., 2012. Combined endocardial and epicardial catheter ablation in arrhythmogenic right ventricular dysplasia incorporating scar dechanneling technique. Circ.: Arrhythmia Electrophysiol. 5 (1), 111–121. Bers, L., 1977. Quasiconformal mappings, with applications to differential equations, function theory and topology. 10.1090/S0 0 02-9904-1977-14390-5 Besl, P., McKay, H., 1992. A method for registration of 3-D shapes. IEEE Trans. Pattern Anal. Mach. Intell. 14 (2), 239–256. Bookstein, F., 1989. Principal warps: thin-plate splines and the decomposition of deformations. IEEE Trans. Pattern Anal. Mach. Intell. 11. doi:10.1109/34.24792. Codreanu, A., Odille, F., Aliot, E., Marie, P., Magnin-Poull, I., Andronache, M., Mandry, D., Djaballah, W., Régent, D., Felblinger, J., de Chillou, C., 2008. Electroanatomic characterization of post-infarct scars comparison with 3-dimensional myocardial scar reconstruction based on magnetic resonance imaging. J. Am. Coll. Cardiol. 52, 839–842. Crum, W., Camara, O., Hill, D., 2006. Generalized overlap measures for evaluation and validation in medical image analysis. IEEE Trans. Med. Imaging 25 (11), 1451–1461. De Craene, M., Tobon-Gomez, C., Butakoff, C., Duchateau, N., Piella, G., Rhode, K.S., Frangi, A.F., 2012. Temporal Diffeomorphic Free Form Deformation (TDFFD) applied to motion and deformation quantification of tagged MRI sequences. In: Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), 7085 LNCS, pp. 68–77. doi:10.1007/978- 3- 642- 28326- 0_7. Duchateau, N., Bijnens, B., D’hooge, J., Sitges, M., 2013. Three-dimensional assessment of cardiac motion and deformation. In: Shiota, T. (Ed.), 3D Echocardiography. CRC Press, pp. 201–213. ˚ Eriksson, A.P., Aström, K., 2012. On the bijectivity of thin-plate splines. In: The Tribute Workshop in Honour of Gunnar Sparr held in Lund, May 8-9, 2008. In: Vol. 6 of Analysis for Science, Engineering and Beyond. Springer, pp. 93–141. doi:10.1007/978- 3- 642- 20236- 0_5. Gu, X., Wang, S., Kim, J., Zeng, Y., Wang, Y., Qin, H., Samaras, D., 2010. Ricci flow for 3D shape analysis. IEEE Patt. Anal. Mach. Intell. 32, 662–677. Hoogendoorn, C., N. Duchateau, N., Sanchez-Quintana, D., Whitmarsh, T., Sukno, F.M., De Craene, M., Lekadir, K., Frangi, A.F., 2013. A high resolution atlas and statistical model of the human heart from multislice ct. IEEE Trans. Med. Imaging 32, 28–44.

J.Fernández-Armenta, Berruezo, A., Andreu, D., Camara, O., Silva, E., Serra, L., Barbarito, V., Carotenutto, L., Evertz, R., Ortiz-Pérez, J.T., Caralt, T.M.D., Perea, R.J., Sitges, M., Mont, L., Frangi, A., Brugada, J., 2013. Three-dimensional architecture of scar and conducting channels based on high resolution ce-cmr: insights for ventricular tachycardia ablation. Circ.: Arrhythmia Electrophysiol. 6.528–237. Karim, R., Ma, Y., Jang, M., Housden, R.J., Williams, S.E., Chen, Z., Ataollahi, A., Althoefer, K., Rinaldi, C.A., Razavi, R., O’Neill, M., Schaeffter, T., Rhode., K.S., 2014. Surface flattening of the human left atrium and proof-of-concept clinical applications. Comput. Med. Imaging Graphics 38 (4), 251–266. Lam, K.C., Gu, X., Lui, L.M., 2014. Genus-one surface registration via Teichmüller extremal mapping. MICCAI 3, 25–32. Li-Fern, H., 2008. Image integration for catheter ablation: searching for the perfect match. Heart Rhythm 5 (4), 536–537. Lui, L.M., Lam, K.C., Yau, S.-T., Gu, X., 2014. Teichmuller mapping (T-map) and its applications to landmark matching registration. SIAM J. Imaging Sci. 391–426. Pinkall, U., Polthier, K., 1993. Computing discrete minimal surfaces and their conjugates. Exp. Math. 2 (1), 15–36. Relan, J., Chinchapatnam, P., Sermesant, M., Rhode, K., Ginks, M., Delingette, H., Rinaldi, C., Razavi, R., Ayache, N., 2011. Coupled personalization of cardiac electrophysiology models for prediction of ischaemic ventricular tachycardia. Interf. Focus.rsfs2010 0 041. Roujol, S., Basha, T., Tan, A., Khanna, V., Chan, R.H., Moghari, M.H., Rayatzadeh, H., Shaw, J.L., Josephson, M.E., Nezafat, R., 2012. Improved multi-modality data fusion of late gadolinium enhancement MRI to left ventricular voltage maps in ventricular tachycardia ablation. IEEE Trans. Biomed. Eng. 60 (5), 1308–1317. Sasaki, T., Miller, C., Hansford, R., Yang, J., Caffo, B., Zviman, M., Henrikson, C., Marine, J., Spragg, D., Cheng, A., Tandr, H., Sinha, S., Kolandaicelu, A., Zimmerman, S., Bluemke, D., Tomaselli, G., Berger, R., Calkins, H., Halperin, H., Nazarian, S., 2012. Myocardial structural associations with local electrograms: a study of post-infarct ventricular tachycardia pathophysiology and magnetic resonance based non-invasive mapping. Circ.: Arrhythmia Electrophysiol. 5, 1081–1090. Smith, N., de Vecchi, A., McCormick, M., Nordsletten, D., Camara, O., Frangi, A., Delingette, H., Sermesant, M., Relan, J., Ayache, N., Krueger, M., Schulze, W., Hose, R., Valverde, I., Beerbaum, P., Staicu, C., Siebes, M., Spaan, J., Hunter, P., Weese, J., Lehmann, H., Chapelle, D., Rezavi, R., 2011. euheart: personalized and integrated cardiac care using patient-specific cardiovascular modelling. Interf. Focus 1, 349–364. Soto-Iglesias, D., Butakoff, C., Andreu, D., Fernández-Armenta, J., Berruezo, A., Sitges, M., Camara, O., 2013. Analyzing electrical patterns in an experimental swine model of dyssynchrony and CRT. In: IEEE Computing in Cardiology Cinference (CinC), vol. 4, pp. 623–626. Tao, Q., Milles, J., Van Huls Van Taxis, C., Lamb, H., Reiber, J.H., Zeppenfeld, K., Van Der Geest, R.J., 2012. Toward magnetic resonance-guided electroanatomical voltage mapping for catheter ablation of scar-related ventricular tachycardia: a comparison of registration methods. J. Cardiovasc. Electrophysiol. 23 (1), 74–80. Tobon-Gomez, C., Zuluaga, M.A., Chubb, H., Williams, S.E., Butakoff, C., Karim, R., Camara, O., Ourselin, S., Rhode, K.S., 2015. Standardised unfold map of the left atrium: regional definition for multimodal image analysis. J. Cardiovasc. Magn. Reson. 17, 41. Trayanova, N.A., Chang, K.C., 2016. How computer simulations of the human heart can improve anti-arrhythmia therapy. J. Physiol. doi:10.1113/JP270532. Tutte, W., 1963. How to draw a graph. Proc. London Math. Soc. 13 (3), 743–768. Ukwatta, E., Arevalo, H., Rajchl, M., White, J., Pashakhanloo, F., Prakosa, A., Herzka, D., McVeigh, E., Lardo, A., Trayanova, N., Vadakkumpadan, F., 2015. Image-based reconstruction of three-dimensional myocardial infarct geometry for patient-specific modeling of cardiac electrophysiology. Med. Phys. 42, 4579–4590. Vaillant, M., Glaunes, J., 2005. Surface matching via currents. In: Information Proccessing in Medical Imaging, series Lecture Notes in Computer Science (LNCS), 3565, pp. 381–392. Vera, S., Ballester, M.A., Gil, D., 2014. Anatomical parametrization for volumetric meshing of the liver. SPIE Med. Imaging 9034. Vera, S., Ballester, M.A., Gil, D., 2015. A novel cochlear reference frame based on the laplace equation. Comput. Assist. Radiol. Surg.: CARS. Yushkevich, P.A., Zhang, H., Gee, J.C., 2006. Continous medical representation for harmonic structures. IEEE Trans. Med. Imaging 25, 1547–1564. Zipes, D.P., Camm, A.J., Borggrefe, M., Buxton, A.E., Chaitman, B., Fromer, M., Gregoratos, G., Klein, G., Moss, A.J., Myerburg, R.J., Priori, S.G., Quinones, M.A., Roden, D.M., Silka, M.J., Tracy, C., 20 06. ACC/AHA/ESC 20 06 guidelines for management of patients with ventricular arrhythmias and the prevention of sudden cardiac death. Europace. 8, 746–837.