Cobiveco: Consistent biventricular coordinates for precise and intuitive description of position in the heart -- with MATLAB implementation
Steffen Schuler, Nicolas Pilia, Danila Potyagaylo, Axel Loewe
CCobiveco: Consistent biventricular coordinates for precise and intuitive description ofposition in the heart – with MATLAB implementation
Ste ff en Schuler a,1, ∗ , Nicolas Pilia a,1 , Danila Potyagaylo b , Axel Loewe a a Institute of Biomedical Engineering, Karlsruhe Institute of Technology (KIT), 76131 Karlsruhe, Germany b EPIQure GmbH, 76139 Karlsruhe, Germany
Abstract
Ventricular coordinates are widely used as a versatile tool for various applications thatbenefit from a description of local position within the heart. However, the practicalusefulness of ventricular coordinates is determined by their ability to meet application-specific requirements. For regression-based estimation of biventricular position, forexample, a consistent definition of coordinate directions in both ventricles is important.For the transfer of data between di ff erent hearts as another use case, the coordinatevalues are required to be consistent across di ff erent geometries. Existing ventricular co-ordinate systems do not meet these requirements. We first compare di ff erent approachesto compute coordinates and then present Cobiveco, a consistent and intuitive biventric-ular coordinate system to overcome these drawbacks. A novel one-way mapping erroris introduced to assess the consistency of the coordinates. Evaluation of mapping andlinearity errors on 36 patient geometries showed a more than 4-fold improvement com-pared to a state-of-the-art method. Finally, we show two application examples under-lining the relevance for cardiac data processing. Cobiveco MATLAB code is availableunder a permissive open-source license.
1. Introduction
Patient and species independent representations of ventric-ular anatomy are a valuable tool for data processing in cardi-ology. Typical applications include the standardized visualiza-tion and regional evaluation of cardiac data, a transfer of databetween di ff erent hearts from di ff erent measurement modali-ties, and the description of local position in the heart. Re-cently, such representations have become particularly importantfor non-invasive localization of the excitation origin using ma-chine learning algorithms (Yang et al., 2018; Zhou et al., 2019).The most popular example of such a representation is the AHAsegmentation from Cerqueira et al. (2002) which divides the left ∗ Corresponding author: Institute of Biomedical Engineering, Karlsruhe In-stitute of Technology (KIT), Fritz-Haber-Weg 1, 76131 Karlsruhe, Germany.Email address: [email protected] (S. Schuler). These authors contributed equally. ventricle (LV) into 17 segments. While easy to apply in prac-tice, it only allows a discrete, coarse-grained representation ofthe LV and does not cover the right ventricle (RV).The approach proposed by Paun et al. (2017) is more generalas it provides a continuous parameterization of both LV andRV. It uses solutions to Laplace’s equation to flatten a ventric-ular bounding surface onto a planar domain and to encode thethickness of anatomical structures on top of this planar domain.Although intended for detailed representation of the ventricu-lar interior (endocardium, trabeculations, papillary muscles), itmay also be applied to the whole myocardial wall. However,this approach does not directly provide an intuitive descriptionof local position and treats LV and RV independently.The universal ventricular coordinates (UVC) introducedby Bayer et al. (2018) o ff er such an intuitive descriptionby defining an apicobasal, a rotational, a transmural and atransventricular coordinate – each of which is defined using so-lutions to Laplace’s equation. UVC thereby o ff er a parameter- a r X i v : . [ ee ss . SP ] F e b S. Schuler et al. (2021) ized description of ventricular position and a similar methodexists for the atria (Roney et al., 2019). Nevertheless, we ar-gue that the UVC system is not consistent in three ways: First,the definition of coordinates in the LV and RV is not symmetric,which causes discontinuities at the junctions of LV and RV. Thiscan be problematic, for instance, for regression-based estima-tion of ventricular position. Second, the definition of boundaryconditions can lead to a non-zero apicobasal coordinate at sin-gularities of the rotational coordinate, which may result in un-defined behavior at these locations. Third, the coordinate valuesthemselves are not very consistent across di ff erent geometries,as solutions to Laplace’s equation are in general not an accu-rate measure of (normalized) distances. These properties of theUVC can lead to errors in mapping data between hearts andother applications that rely on a consistent description of localposition. This especially pertains to the validation of electrocar-diographic imaging (Cluitmans et al., 2018). Here, the transferof potentials or activation times from a geometry obtained us-ing intracardiac mapping onto a tomography-derived geometry,which is used for inverse reconstructions, is usually needed. In-consistencies between coordinates computed on these two ge-ometries can cause problematic artifacts in the mapped signals.In this work, we propose a new coordinate system for biven-tricular geometries that builds upon the ideas from Bayer et al.(2018) and previous works but removes inconsistencies andtherefore reduces mapping errors. We start by defining desir-able properties for such a coordinate system. Then, we explainthe underlying concept for consistent coordinate directions andcompare di ff erent approaches for computing the actual coordi-nate values based on partial di ff erential equations (PDE). Hav-ing identified a suitable approach, we provide a detailed de-scription of the new coordinate system, called Cobiveco . Fi-nally, we compare Cobiveco with UVC by evaluating mappingand linearity errors and present application examples.
2. Methods
Based on the use cases mentioned in the introduction, thefollowing properties are considered desirable for a biventricularcoordinate system: • Bijective:
Each coordinate tuple corresponds to exactlyone point in the heart. • Continuous:
Coordinates have no jumps. • Normalized:
Coordinates range between 0 and 1. • Complete:
Each tuple in this range represents a valid posi-tion, i.e., the coordinate space has no holes. • Linear:
Coordinates change linearly in space, i.e., thegeodesic distance traveled when changing one coordinate,while keeping all others fixed, is proportional to the changein this coordinate. • Consistent parameterization:
The underlying parameteri-zation is the same for both ventricles. • Consistent landmarks:
Clear anatomical landmarks arerepresented by the same coordinates across di ff erenthearts. In particular, landmarks used to construct the coor-dinate system are robust to variations in shape.Note that normalized and linear coordinates can in general(for arbitrary shapes) not also be orthogonal. The resulting co-ordinate system will not preserve angles, but it will preservedistances in each of the coordinate directions. Following the UVC approach from Bayer et al. (2018), ourchoice of coordinate directions is inspired by prolate spheroidalcoordinates as used in Costa et al. (1996) to parameterize anidealized LV geometry using an apicobasal, a rotational anda transmural coordinate. The goal of UVC and Cobiveco isto find a “generalization” for biventricular geometries of arbi-trary shape. To this end, one more transventricular coordinateis needed that distinguishes between LV and RV.The left panel of Fig. 1 illustrates the basic concept for thesefour coordinates within the UVC system. Here, the transven-tricular boundary is chosen such that the entire septum belongsto the LV. While this choice is anatomically intuitive and mightbe most useful for applications focusing on the LV, it leads toundesired properties of the coordinates: The transmural and therotational coordinates are discontinuous at the transventricularboundary and the ranges of the rotational coordinate are di ff er-ent in the LV and the RV ( − π to π vs. − π/ π/ UVC Cobiveco −11 0 001 10 0; 1 01 2/3Transventricular coordinateTransmural coordinateRotational coordinateApicobasal coordinate Transventricular boundaryImaginary axis of rotationRotational singularityRotational boundaryand apicobasal origin± π π /2 π /2.5− π /2− π /2.5 01 01 1 Fig. 1. Di ff erent underlying concepts for coordinate directions and boundaryvalues within the UVC ( left ) and the suggested coordinate system Cobiveco( right ). A basal cross-section in long-axis direction and a central cross-sectionin anterior-posterior direction are shown.. Schuler et al. (2021) 3 Solving PDEs can be an e ffi cient and elegant way to computecoordinate values. However, the type of PDE and the boundaryconditions need to be chosen with care. In this section, dif-ferent PDE-based approaches to compute coordinate values arepresented and compared to each other to choose the most ade-quate one. Solutions to Laplace’s equation (in the following just“Laplace solutions”) in between two boundaries with Dirichletconditions are one obvious approach and were utilized in Bayeret al. (2018). Nevertheless, they are not necessarily a goodchoice as their linearity severely depends on the width of thedomain between these boundaries. This follows directly fromthe divergence theorem: As the Laplace equation requires thedivergence of the gradient to be zero and the flux of the gra-dient field through lateral parts of the outer surface is alreadyzero due to zero Neumann boundary conditions, the (signed)fluxes of the gradient field through any two cross-sectional sur-faces have to compensate each other. This implies that smallergradients occur in wider regions and vice versa, which makesLaplace solutions an unreliable measure of normalized distancebetween boundaries. Therefore, they should not directly beused to define coordinate values. In the UVC method, this issuegets obvious for the apicobasal Laplace solution in between asmall apical and a large basal boundary, where it led the authorsto normalize the resulting apicobasal coordinate using the val-ues on the shortest geodesic path between apex and base (Bayeret al., 2018). However, it can also lead to substantial distortionsof the transmural and rotational coordinates.To demonstrate the e ff ect on the rotational coordinate, wecreated two ellipsoidal geometries resembling the LV free wall:One with uniform and one with non-uniform wall thickness asexpected in reality (Fig. 2, top). The Laplace solution u be-tween the two boundary surfaces S and S is depicted in thefirst column of Fig. 2: ∆ u = u ( S ) = u ( S ) = g and g with respect to asingle boundary can be obtained: (cid:107)∇ g (cid:107) = g ( S ) = (cid:107)∇ g (cid:107) = g ( S ) = g = g g + g (4)However, the result shows a very inhomogeneous distributionof contour lines, even for the case of uniform wall thickness Laplace Eikonal Trajectory u g = g g + g d = d d + d g g d d Poisson p = p p + p p p Uniform wall thickness
Non - uniform wall thickness S S N o n - un i f o r m U n i f o r m w a llt h i c kn e ss … Fig. 2.
Comparison of four approaches to compute a rotational coordinate.
Top:
Ellipsoidal geometries with boundary surfaces used as input.
Bottom:
Resultsfor the di ff erent approaches. For the Eikonal approach, several cusps occur atthe lateral surfaces. The green arrow marks one of them. (second column of Fig. 2). Furthermore, the contour lines oftenhave cusps (green arrow). The reason is that g and g representdistances along di ff erent, non-bijective trajectories between S and S , which makes the normalization according to (4) invalid.Another strategy to reduce the non-linearities of the Laplacesolution is to compute its gradient field, normalize it to unitlength and “integrate it back” by solving Poisson’s equation: t = ∇ u (cid:107)∇ u (cid:107) (5) ∆ p = ∇ · t with p ( S ) = ∆ p = −∇ · t with p ( S ) = p for the uniform but not for the non-uniform case (third column of Fig. 2). The problem here is thatalthough the trajectories along t are bijective, the trajectoriesalong ∇ p and ∇ p are not anymore, due to the divergence op-erator in (6), (7). Therefore, Poisson’s equation is not an appro-priate way to integrate t for our purpose. A further way is bysolving the “trajectory distance equation” originally proposedfor obtaining a symmetric measure of tissue thickness by Yezziand Prince (2003): ∇ d · t = d ( S ) = −∇ d · t = d ( S ) = S. Schuler et al. (2021)
These systems of linear equations are overdetermined, becausethere are more elements than nodes in a tetrahedral mesh.Hence, they are solved in a least-squares sense. As the gra-dient fields of the trajectory distances d and d themselvesnow match t and − t , respectively (and not only their diver-gences), they allow for normalization as in (4). The normalizedresult d is shown in the last column of Fig. 2 and exhibits thedesired behavior even for the non-uniform case.We conclude that normalized distances obtained by solving thetrajectory distance equation are well suited to define coordinatevalues and should be preferred over Laplace’s, Poisson’s or theEikonal equation. In each of the following subsections, we will describe one ofthe eight steps involved in the computation of coordinate valuesaccording to Cobiveco. We use V for denoting volumes, S forsurfaces, C for curves and x for points. u is used for Laplacesolutions and d for (relative) trajectory distances. The mathe-matical description is paralleled by a purely verbal description.To increase readability, we use verbal terms introduced in italicfont instead of the corresponding mathematical symbols when-ever possible in the text. The computational process on discretemeshes is described in the text and using pictures, while themathematical notation refers to the continuous case.We provide an open-source MATLAB implementation ofCobiveco for tetrahedral meshes ( https://github.com/KIT-IBT/Cobiveco ) under the Apache License 2.0. For moredetails about the implementation, the reader is referred to thecode. To solve partial di ff erential equations, we use the discreteLaplace and gradient operators from the gptoolbox (Jacobson,2018), which are equivalent to first order finite element dis-cretization (Jacobson, 2013). For general geometry processing(thresholding, surface extraction, connectivity filtering, isocon-tour computation, etc.), the VTK library (Schroeder et al., 2006)is used, for which we have developed and provide a MEX in-terface called vtkToolbox ( https://github.com/KIT-IBT/vtkToolbox ). Implicit domain remeshing (isovalue discretiza-tion) is performed using mmg3d (Dapogny et al., 2014). Cobiveco requires a biventricular volume V with exactly oneorifice at the base of each ventricle. If there are bridges betweenthe tricuspid valve and the RV outflow tract or between the mi-tral valve and the LV outflow tract, they have to be removed.To yield consistent results across di ff erent geometries, the baseof the heart should be truncated at comparable heights. Apartfrom the volume mesh, four boundary surfaces as depicted inFig. 3 are needed as input: a basal surface S Base , an epicardialsurface S
Epi , an
LV endocardial surface S LV , and an RV endo-cardial surface S RV .We provide utilities for semi-automatic clipping at the base, re-moval of bridges and extraction of these boundary surfaces aspart of the Cobiveco code. S Base S Epi S LV S RV Fig. 3.
Boundary surfaces required as input: Basal surface, epicardial surface,LV endocardial surface and RV endocardial surface.
To compute the transventricular coordinate, we first solveLaplace’s equation with boundary conditions of 0 at the RV en-docardium and 1 at the LV endocardium (Fig. 4, left): ∆ u v ( V ) = u v ( S RV ) = u v ( S LV ) = transven-tricular coordinate v with binary values (Fig. 4, right): v = round( u v ) (11) u v v Fig. 4.
Computation of the transventricular coordinate.
Left:
Laplace solution.
Right:
Final coordinate. The geometry was clipped for visualization.
For being able to apply boundary conditions exactly at theboundary between the LV and the RV, we perform isovalue dis-cretization at u v = .
5, which yields mesh 1 (Fig. 5). Thismeans that the original tetrahedral mesh is remeshed, such thatthere are nodes directly on the boundary between the two ven-tricles. u v < 0.5 > 0.50.5 Original mesh Mesh 1 Fig. 5.
Close-up of the clipped mesh at the anterior interventricular junctionbefore ( left ) and after ( right ) isovalue discretization at u v = . From mesh 1, we extract all tetrahedron faces composing thisboundary, which results in a septal surface S
Sept (Fig. 6, left): S Sept = (cid:8) x ∈ V | u v ( x ) = . (cid:9) (12) . Schuler et al. (2021) 5 Similarly, we can extract a septal curve C
Sept from the corre-sponding epicardial surface (Fig. 6, right): C Sept = (cid:8) x ∈ S Epi | u v ( x ) = . (cid:9) (13) S Sept C Sept
Fig.
6. Septal surface ( left ) and septal curve ( right ). To obtain the transmural coordinate, we first compute aLaplace solution that is 0 at the epicardial and the septal sur-face and 1 at the LV and RV endocardial surfaces (Fig. 7, left): ∆ u m ( V ) = u m ( S Epi ∪ S Sept ) = u m ( S LV ∪ S RV ) = d m along the gradient ofthis Laplace solution in both directions, i.e., starting from theepicardium and starting from the endocardium: ∇ d m , Epi · t m = d m , Epi ( S Epi ∪ S Sept ) = −∇ d m , Endo · t m = d m , Endo ( S LV ∪ S RV ) = t m = ∇ u m (cid:107)∇ u m (cid:107) (17)The relative trajectory distance with respect to the epicardiumis then defined as the transmural coordinate m (Fig. 7, right): m = d m , Epi d m , Epi + d m , Endo (18)Equations (14)-(18) are solved on mesh 1. Linear interpolationis used to transfer the transmural coordinate back to the originaltetrahedral mesh. m u m Fig. 7.
Computation of the transmural coordinate.
Left:
Laplace solution.
Right:
Final coordinate. The geometry was clipped for visualization.
For the rotational and apicobasal coordinate, a consistent androbust definition of an epicardial apex point is essential. Asthis point will be used to define the apex for both ventricles, itshould lie at the center between the two ventricles. Therefore,possible points are restricted to the septal curve in Fig. 6 (right).The most straightforward choice would be the point on this sep-tal curve with the maximum distance to the basal surface. How-ever, this definition would not be very robust, as the positionalong the septal curve would largely depend on its local shapeand smoothness in the apex region. To yield an intuitive courseof the rotational coordinate in the LV, the apex point shouldfurthermore be centered with the LV in anterior-posterior direc-tion. For this reason, we decided to take a more global approachthat relies on the definition of orthogonal heart axes as depictedin the left half of Fig. 8.The long axis v LongAx is defined as the unit vector “most or-thogonal” to the normals of the LV endocardial surface, as mea-sured by the dot product with all triangle normals n LV : v LongAx = arg min v ∈ R (cid:107) v · n LV (cid:107) p (19)Here, (cid:107) · (cid:107) p denotes the p -norm across all surface triangles. Wechose p = . p = p =
2. Similar valueswork as well. Problem (19) is solved using the Nelder-Meadalgorithm. To assure that the long axis is directed from basetowards apex, its dot product with the vector pointing from thecentroid of the basal surface to the centroid of the LV endocar-dial surface is evaluated and the long axis is flipped accordingly.The definition of the left-right axis is based on fitting a planeto the septal surface in Fig. 6 (left). As the septal surface maybecome strongly curved near the interventricular junctions, par-ticularly at the anterior side, a two-step process is used to onlytake into account the central part of the septal surface.In the first step, principal component analysis is applied to thepoints on the entire septal surface. The third principal compo-nent represents the normal vector of the best-fitting plane andis defined as v LR , Entire . Here, the vector pointing from the cen-troid of the LV endocardial surface to the centroid of the RVendocardial surface is used as reference to assure that this vec-tor is directed from left to right. The distance in the directionof v AP = v LongAx × v LR , Entire is then used to truncate the septalsurface by 20 % and 10 % at the anterior and posterior side, re-spectively, which yields the truncated septal surface S
SeptTrunc : S SeptTrunc = (cid:8) x ∈ S Sept | x · v AP > P ( x · v AP ) and (20) x · v AP < P ( x · v AP ) (cid:9) (21)Here, P q denotes the q th percentile.In the second step, the final left-right axis v LeftRightAx is obtainedby computing the third principal component v LR , Trunc of pointson the truncated septal surface and orthogonalizing it with re-spect to the long axis: v LeftRightAx = v LR , Trunc − ( v LR , Trunc · v LongAx ) v LongAx (22)The anterior-posterior axis v AntPostAx is finally defined as: v AntPostAx = v LongAx × v LeftRightAx (23)
S. Schuler et al. (2021) v AntPostAx C SeptAnt C SeptPost v LongAx v LeftRightAx S SeptTrunc x LV x Apex x Center
Fig. 8.
Left:
Truncated septal surface (dark yellow) and heart axes (blue, redand green arrows).
Right:
Steps to locate the apex point: The LV centroid (reddot) is projected (red line) onto the plane defined by the left-right axis and theseptal centroid, which yields a global center point (blue dot). A line in long axisdirection (blue) is starting at this center point and the point of the septal curveclosest to this line is identified as apex point (yellow dot). The apex point splitsthe septal curve into an anterior part (cyan) and a posterior part (magenta).
The three heart axes are then used to find the apex point.This is illustrated in the right half of Fig. 8. First, a global center point x Center is obtained by projecting the centroid x LV ofthe LV endocardial surface onto the plane perpendicular to theleft-right axis that passes through the centroid x SeptTrunc of thetruncated septal surface (red line): x Center = x LV + (( x SeptTrunc − x LV ) · v LeftRightAx ) v LeftRightAx (24)The apex point x Apex is then located as the point of the septalcurve with the smallest distance to the line in long axis directionstarting at this center point (blue line): x Apex = arg min x ∈ C Sept (cid:107) x + r LongAx v LongAx − x Center (cid:107) (25)with r LongAx = ( x Center − x ) · v LongAx > anteriorseptal curve C SeptAnt and a posterior septal curve C
SeptPost . Computing a rotational coordinate by solving a PDE requiresto define at least two surfaces for assigning boundary con-ditions. For consistency, these surfaces should be based onanatomical landmarks that can be identified reliably on di ff erentgeometries. As we furthermore aim for a rotational coordinatethat is symmetric in both ventricles and that allows to distin-guish between the septum and the free walls, the anterior andposterior junctions between the septum and both free walls are anatural choice for such landmarks. To obtain boundary surfacesrepresenting these two junctions, we first compute a Laplacesolution that is 0 on the epicardial surface and 1 on the septal u Ridge u Ridge < 0.5 > 0.50.5 Mesh 1 Mesh 2
Fig. 9.
Upper half: “Ridge” Laplace solution. The geometry was clipped forvisualization.
Lower half:
Close-up of the mesh at the anterior septal junctionbefore ( left ) and after ( right ) isovalue discretization at u Ridge = . surface (see upper half of Fig. 9): ∆ u Ridge ( V ) = u Ridge ( S Epi \ S Sept ) = u Ridge ( S Sept ) = u Ridge = .
5, whichyields mesh 2 (lower half of Fig. 9). Note that in the boundaryconditions of (26), the epicardial points of the septal surfaceare excluded from the epicardial surface to obtain disjoint meshregions (no common nodes) for the left and right free walls afterremeshing. From mesh 2, we extract the volume V Free coveringboth free walls, the volume V Sept covering the septal wall, and a ridge surface S
Ridge : V Free = (cid:8) x ∈ V | u Ridge ( x ) ≤ . (cid:9) (27) V Sept = (cid:8) x ∈ V | u Ridge ( x ) ≥ . (cid:9) (28) S Ridge = (cid:8) x ∈ V | u Ridge ( x ) = . (cid:9) (29)The shortest Euclidean distance to the anterior and posteriorseptal curves is used to split the ridge surface into an anteriorridge surface S RidgeAnt and a posterior ridge surface S
RidgePost in a nearest neighbor manner. Furthermore, a transmural apexcurve C
Apex is obtained between these two surfaces (Fig. 10): S RidgeAnt = (cid:8) x ∈ S Ridge | r Ant ( x ) < r Post ( x ) (cid:9) (30) S RidgePost = (cid:8) x ∈ S Ridge | r Ant ( x ) > r Post ( x ) (cid:9) (31) C Apex = (cid:8) x ∈ S Ridge | r Ant ( x ) = r Post ( x ) (cid:9) (32)with r Ant ( x ) = min y ∈ C SeptAnt (cid:107) x − y (cid:107) , r Post ( x ) = min y ∈ C SeptPost (cid:107) x − y (cid:107) As there are no nodes that exactly fulfill r Ant ( x ) = r Post ( x ), theclosest nodes on the anterior ridge surface define the apex curvein the discrete mesh. . Schuler et al. (2021) 7 S RidgeAnt S RidgePost C Apex V Free V Sept
Fig. 10.
Anterior and posterior ridge surfaces and apex curve. The apex curveruns from the epicardium to the LV and RV endocardium.
The relative trajectory distance between the posterior and an-terior ridge surfaces is used to define the rotational coordinate(Fig. 11, top-left). It is computed separately within the freewalls and the septum: d r ( V Free ) = d r , Post ( V Free ) d r , Post ( V Free ) + d r , Ant ( V Free ) (33) d r ( V Sept ) = d r , Post ( V Sept ) d r , Post ( V Sept ) + d r , Ant ( V Sept ) (34)where d r , Post and d r , Ant are given by: ∇ d r , Post ( V Free ) · t r ( V Free ) = d r , Post ( S RidgePost ) = −∇ d r , Ant ( V Free ) · t r ( V Free ) = d r , Ant ( S RidgeAnt ) = ∇ d r , Post ( V Sept ) · t r ( V Sept ) = d r , Post ( S RidgePost ) = −∇ d r , Ant ( V Sept ) · t r ( V Sept ) = d r , Ant ( S RidgeAnt ) = t r is not based on the gradient of a rota-tional Laplace solution but on the cross product of the gradientsof the transmural coordinate m (cid:48) and an apicobasal Laplace so-lution u a : t r = ∇ m (cid:48) (cid:107)∇ m (cid:48) (cid:107) × ∇ u a (cid:107)∇ u a (cid:107) (39)The transmural coordinate is inverted in the RV to get coherentgradients in the septum and opposite directions of rotation inboth ventricles: m (cid:48) ( x ) = m ( x ) , v ( x ) = − m ( x ) , v ( x ) = m (cid:48) from mesh 1 to mesh 2. The apicobasalLaplace solution (Fig. 11, top-right) is computed directly onmesh 2 and is 0 at the apex curve and 1 at the basal surface: ∆ u a ( V ) = u a ( C Apex ) = u a ( S Base ) = rotational coordinate r (Fig. 11, bottom) is obtainedby flipping, scaling and shifting the relative trajectory distances: r ( V Free ) = d r ( V Free ) (42) r ( V Sept ) = + (cid:0) − d r ( V Sept ) (cid:1) (43)Based on average geometrical proportions and in accordancewith the ratio of two septal and four free wall segments in theAHA scheme (Cerqueira et al., 2002), the scaling factors werechosen such that the septum covers one third and the free wallstwo thirds of the total range [0 , / rotational sinecoordinate r sin and a rotational cosine coordinate r cos : r sin = sin(2 π r ) (44) r cos = cos(2 π r ) (45)This trick is used for linear interpolation back onto the originalmesh, where the following inverse transform is applied: r = atan2( r sin , r cos ) / (2 π ) , r sin ≥ r sin , r cos ) / (2 π ) + , r sin < r d r r u a Fig. 11.
Computation of the rotational coordinate.
Top-left:
Relative trajectorydistance.
Top-right:
Apicobasal Laplace solution.
Bottom:
Final coordinate.The geometry shown on the left was clipped for visualization. S. Schuler et al. (2021)
Although trajectory distances between two boundary sur-faces are used to define the transmural and the rotational co-ordinate, this approach is not well suited for the apicobasal co-ordinate. The reason is that in this case, we are looking fora normalized distance between the two-dimensional basal sur-face and the one-dimensional apex curve. Due to the di ff erentdimensionality of boundaries, trajectories starting at di ff erentpoints on the basal surface may end at the same point on theapex curve, which leads to contradicting values of the trajectorydistance. Therefore, a di ff erent approach is used: Apicobasalcurves are obtained by extracting isocontours at discrete valuesof the transmural and the rotational coordinate and the normal-ized distance along these curves is determined.We start by extracting isosurfaces of the transmural coordinatefrom mesh 1. The following 10 isovalues are used to obtain anequidistant sampling: m ∈ (cid:110) , , , . . . , (cid:111) (47)This results in 20 disjoint isosurfaces S i ( i = , , . . . , ffi ciently fine sampling that captures the septal junctions at r = / r = r ∈ (cid:110) , , , . . . , (cid:111) (48)This results in 1920 isocurves C i , j ( i as in S i , j = , , . . . , S i and one in-dividual apex point x i is determined for each S i by finding theminimum of the Laplace solution. Then, the curves are trun-cated by excluding points within a radius ε to the respectiveapex point: C Trunc i , j = (cid:8) x ∈ C i , j | (cid:107) x − x i (cid:107) > ε (cid:9) with x i = arg min x ∈ S i u a ( x ) (49)An ε of three times the mean edge length of the original meshwas found to be su ffi cient to ensure disjoint curves.To obtain smooth curves with well-defined apical start points,the corresponding x i is re-added to each C Trunc i , j and a cubicsmoothing spline fit (Reinsch, 1967) is used to resample eachcurve at 100 equidistant nodes along the curve. This yieldsthe spline curves C Spline i , j . The extent of smoothing is deter-mined such that the root-mean-square deviation (RMSD) fromthe original points equals 0 . x i , a 100-fold weight is used for this point. The normalized distance a Spline along each spline curve is then computed as therelative cumulative sum of Euclidean distances between neigh-boring nodes on this curve, starting at x i . The result can be seenat the bottom-left of Fig. 12.Finally, Laplacian extrapolation is used to obtain the apicobasalcoordinate a on the original mesh from a Spline : a = arg min a (cid:16) (cid:107) Ra − a Spline (cid:107) + λ (cid:107) La (cid:107) + η (cid:107) Ea − (cid:107) (cid:17) (50) = ( R T R + λ L T L + η E T E ) − ( R T a Spline + η E T )Here, the vector a Spline contains a Spline at all nodes of the splinecurves and a contains a at all nodes of the volume mesh. R isa matrix that linearly interpolates from the nodes of the volumemesh onto the nodes of the spline curves and L is the Laplacianoperator of the volume mesh. The smoothing parameter λ isdetermined using the secant method, such that the RMSD be-tween Ra and a (cid:48) equals 25 % of the mean edge length of thevolume mesh divided by the mean length of the curves. Thelast term in (50) forces the extrapolated values to 1 at the base. E extracts the values at the basal surface of the volume meshand η is chosen to yield an equal weighting with the first term: η = (cid:16) number of nodes on the spline curvesnumber of nodes on the basal surface (cid:17) (51)The final apicobasal coordinate after extrapolation is depictedat the bottom-right of Fig. 12. a a Spline C Spline i , j C Trunc i , j u a S i on Fig. 12.
Computation of the apicobasal coordinate.
Upper half:
ApicobasalLaplace solution on the LV and RV isosurfaces for one out of 10 transmuralvalues and corresponding isocurves for all 96 rotational values (black lines) af-ter truncation and spline fitting.
Lower half:
Normalized distance on all 2 · · left ) and final coordinate on the original mesh ( right ).. Schuler et al. (2021) 9 To transfer scalar data using Cobiveco, we construct a ma-trix M A ← B that maps from the nodes of a source mesh B to thenodes of a target mesh A (or a target point cloud). The princi-pal mapping procedure is similar to the one described in Bayeret al. (2018) – with three advancements: • The discontinuity of the rotational coordinate is com-pletely avoided by transforming it into the continuous sineand cosine coordinates using (44) and (45). • Mesh dependent instead of fixed scaling factors are usedto yield scaled coordinates that show a comparable changeper unit length in Euclidean space. To this end, the max-imum coordinate di ff erence of m , r sin , r cos and a betweenany two nodes of each tetrahedron is computed for B . Thecoordinates in both A and B are then divided by the me-dian value of the respective maximum di ff erences. As thisis not possible for the binary coordinate, v is instead mul-tiplied by the bounding box diagonal and divided by themean edge length of B . • The rotational coordinates are additionally scaled as afunction of the apicobasal coordinate. This is important toassure a well-defined mapping at the rotational singulari-ties and to account for the decreasing circumference of theventricles towards the apex. As the rotational coordinatebecomes undefined at the apex curve (see Fig. 10 and 11),its weighting should become zero for a =
0. As the cir-cumference of the ventricles is roughly proportional to thesquare root of the apicobasal coordinate, √ a is chosen asscaling function for r sin and r cos .The mapping functionality is also implemented in MATLAB.The user can choose between a mapping using linear or nearest-neighbor interpolation. For linear interpolation, the computa-tion of M A ← B consists of five steps:1. The coordinates are scaled as described above.2. For each point in A , the tetrahedron centroid in B with theclosest ventricular coordinates is found. This is done usinga nearest neighbor search with a k-d tree and an Euclideandistance metric.3. For each centroid found in step 2, all centroids within apredefined search radius are found. This is done separatelyfor the left and right ventricle using a range search with ak-d tree and an Euclidean distance metric. A search radiusof two mean edge lengths of B was found to be su ffi cientand used in this work.4. For each point in A , we iterate over the tetrahedrons cor-responding to the respective centroids found in step 3. Foreach tetrahedron, we compute the barycentric coordinatesthat reproduce the ventricular coordinates of A . The tetra-hedron with the smallest maximum absolute deviation ofbarycentric coordinates from 0 . M A ← B . For nearest-neighbor interpolation, steps 2-4 are replaced by di-rectly finding nodes instead of tetrahedron centroids and M A ← B is made up of ones instead of barycentric coordinates (step 5).
3. Evaluation
Two sets of biventricular geometries were used to evaluateCobiveco: Geometries created using a statistical shape model(SSM) and imaged patient geometries.
The mean shape of the SSM from Bai et al. (2015); de Mar-vao et al. (2014) was used as a representative geometry to eval-uate mapping errors and 1000 quasi-random instances of thismodel were used to assess the computational robustness of Co-biveco. This SSM was created from more than 1000 magneticresonance images. Originally, it consists of disconnected sur-faces of the LV endo- and epicardium and the RV blood pool.To derive a model that can be used to compute coordinates, weextruded the RV blood pool by 3 mm to obtain an RV epicardialsurface and merged all surfaces to form one closed surface ofthe biventricular myocardium. This surface was tetrahedralizedand the 100 principal components and variances were interpo-lated to the nodes of the volume mesh. The adapted model isavailable at https://doi.org/10.5281/zenodo.4419784 .Mesh statistics for the mean shape depicted in Fig. 3 can befound in the first row of Table 3. The 1000 quasi-random in-stances were created by drawing the weights of the 100 shapemodes from a uniform distribution within bounds of ±
36 patient geometries were used for a comparison of Co-biveco and UVC under realistic conditions. These geometrieswere acquired as part of validation studies (Revishvili et al.,2015; Chmelevsky et al., 2018) for electrocardiographic imag-ing (ECGI), which adhered to the Declaration of Helsinki andwere approved by the Institutional Review Board of AlmazovNational Medical Research Center in Saint Petersburg, Russia.Written informed consent was obtained from each patient. Car-diac computed tomography (CT) images were obtained frompatients with implanted pacemakers and segmented in a semi-automatic manner with the software of the Amycard 01C EPsystem (EP Solutions SA, Yverdon-les-Bains, Switzerland). Asthis system uses relatively coarse triangle meshes suitable forECGI (edge lengths of 5 to 10 mm), they were first remeshedwith Instant Meshes (Jakob et al., 2015) and then tetrahedral-ized with Gmsh (Geuzaine and Remacle, 2009). Some geome-tries included large parts of the aorta and the pulmonary artery.To yield consistent inputs for the computation of coordinates,we clipped all meshes at the base (where the LV outflow tractintersects the septal plane) and removed the bridge at the baseof the RV. All 36 geometries are shown in Fig. 22 and meshstatistics are given in Table 3. et al. (2021)
For a comparison of Cobiveco with UVC, we also computedUVC coordinates for the mean shape of the SSM and all patientgeometries. To this end, we reimplemented the UVC method inMATLAB according to the description in Bayer et al. (2018).This implementation is also accessible at https://github.com/KIT-IBT/Cobiveco . The UVC method was providedwith the most comparable inputs: The epicardial apex pointidentified by Cobiveco was used as “user-defined” apex pointand the part of the RV endocardial surface with r ∈ [2 / , ν, ρ, φ, z ) were transformed into coordinates ( v (cid:48) , m (cid:48) , r (cid:48) , a (cid:48) ) thatcover the same ranges as the corresponding Cobiveco coordi-nates (see Fig. 1): v (cid:48) = + ν (52) m (cid:48) = − ρ (53) r (cid:48) = + π atan2 (cos φ, sin φ ) , ν = − ∧ | φ | > π/ + π atan2 (cos φ, sin φ ) , ν = − ∧ | φ | ≤ π/ + π φ, ν = a (cid:48) = z (55) One quantitative way to evaluate the ventricular coordinatesis to use them to map the Euclidean coordinates of one heartto a second heart and then back again to the first heart. A de-viation between the original and the mapped Euclidean coordi-nates can then be computed on the first heart. This “two-wayerror” was used in Bayer et al. (2018). However, it only re-flects errors due to non-bijectivity and interpolation. To cap-ture errors due to inconsistencies in the ventricular coordinatevalues across di ff erent geometries, one has to map only in onedirection, i.e., from the first heart to the second heart, and thencompute the deviation to the ground truth on the second heart.However, no real ground truth is available, because no error-free reference method exists to determine the anatomical pointcorrespondences between both hearts. To overcome this prob-lem, we used the ventricular coordinates themselves to createa synthetic “mean heart geometry” for which the ground truthis known by construction. The resulting “one-way error” there-fore reflects the self-consistency of the ventricular coordinates.We first derive mathematical expressions for the two- and one-way errors and then illustrate the one-way error on a simpletwo-dimensional example, before we explain the evaluation onthe actual test geometries.Let k ∈ N N × denote the node indices, X ∈ R N × the Eu-clidean coordinates and V ∈ R N × the ventricular coordinatesof a heart mesh, with N being the number of nodes. We intro-duce a function χ A that maps from node indices to the Euclideancoordinates of a heart A : X A = χ A ( k A ) (56)The relation between Euclidean and ventricular coordinates ofa heart A is represented by a function ϕ A : V A = ϕ A ( X A ) (57) Using this function, the process of mapping Euclidean coordi-nates from a heart B to a heart A can be expressed as: Φ A ← B ( X A ) = ϕ − B ( ϕ A ( X A )) (58)The right side of this equation has to be interpreted in the fol-lowing way: As we want to get values at the nodes of A , weplug X A into ϕ A to obtain the ventricular coordinates V A . Thenwe extract the Euclidean coordinates of B at points with thesame ventricular coordinates by applying ϕ − B to V A . In prac-tice, this process is implemented using the mapping matrix fromsection 2.5 (with linear interpolation): Φ A ← B ( X A ) = M A ← B X B (59)Having introduced these notations, the two-way error for themapping sequence “ A to B and back to A ” can be written as: e two-way AB = (cid:107) X A − (cid:101) X A (cid:107) col with (cid:101) X A = Φ A ← B ( Φ B ← A ( X B )) (60)Here, (cid:107) · (cid:107) col denotes the 2-norm along the column dimension.The novel one-way error between A and B with respect to A isdefined as: e one-way AB = (cid:107) X C − (cid:101) X C (cid:107) col (61)with X C = (cid:0) X A + Φ A ← B ( X A ) (cid:1) (62) (cid:101) X C = Φ C ← A ( χ A ( k C )) = ϕ − C ( ϕ A ( χ A ( k C ))) (63) X C is the average of the Euclidean coordinates at the nodes of A and the Euclidean coordinates at the corresponding points in B . Together with the mesh connectivity of A , it forms the meanheart geometry C . As the node indices of A and C are identical,we can copy the ventricular coordinates computed on A to C .This yields the “ground truth” coordinates V truth C represented by ϕ A ( χ A ( k C )) in (63). A new set of coordinates “to be evaluated” V C is then computed on C and Euclidean coordinates (cid:101) X C are ex-tracted at points where these ventricular coordinates equal theground truth coordinates. This is represented by ϕ − C ( V truth C ). Ifthe coordinate system is consistent across di ff erent geometries,the norm of the di ff erence in (61) should be as small as possi-ble. As the mapping between A and C covers only half the waybetween A and B , the norm in (61) is multiplied by two.Fig. 13 illustrates the principle of the one-way error usinga simple two-dimensional example. Here, the hearts A and B are represented by a star-shaped and a flower-shaped geome-try, respectively. For both geometries, “ventricular” coordinatesare computed analogously to the definition of the apicobasaland the rotational coordinates in Cobiveco and UVC (black ar-rows). In this example, averaging both geometries using theCobiveco coordinates yields an almost circular mean geometry,while the UVC coordinates lead to a spiky mean geometry (bluearrows). Two sets of ventricular coordinates are then obtainedon the respective mean geometry: The coordinates to be evalu-ated are computed independently (red arrow), while the groundtruth coordinates are copied from A (green arrows). For Co-biveco, these two sets of ventricular coordinates look very sim-ilar, while larger di ff erences can be seen for UVC, especially for . Schuler et al. (2021) 11 V B = φ B ( X B ) C ob i v e c o U V C V A = φ A ( X A ) C ob i v e c o U V C V C = φ C ( X C ) C ob i v e c o U V C Coordinates to be evaluated C ob i v e c o U V C Ground truth coordinates V truth C = φ A ( χ A ( k C )) X A Target A X B Source B X C = ( X A + Φ A ← B ( X A ) ) C ob i v e c o U V C Mean C ˜ X C = φ −1 C ( V truth C ) C ob i v e c o U V C a , r e y AB = ∥ X C − ˜ X C ∥ c o l on e - w a Fig. 13.
Two-dimensional example illustrating the computation of the one-way error in (61)-(63). We start with the Euclidean coordinates X B and X A of a sourcegeometry B and a target geometry A , for which we compute the ventricular coordinates V B and V A (black arrows). These ventricular coordinates are used to map X B to X A and the result is averaged with X A , yielding a mean geometry C represented by X C (blue arrows). For this geometry, ventricular coordinates V C can becalculated (red arrow). Furthermore, we can copy V A to X C to obtain V truth C (green arrows). The mapping between V C and (cid:101) X C is then used to transfer V truth C (yellowarrows) into Euclidean space, which yields (cid:101) X C (yellow arrows). Finally, X C is compared with (cid:101) X C (gray arrows). the apicobasal coordinate. To quantify the Euclidean error as-sociated with this inconsistency of ventricular coordinates, weuse the mapping between the coordinates to be evaluated andthe Euclidean coordinates of the mean geometry to transfer theground truth coordinates into Euclidean space (yellow arrows).The result is then compared with the Euclidean coordinates ofthe mean geometry (gray arrows).For the actual test geometries, two- and one-way errors arealways computed between the mean shape of the SSM and oneof the patient geometries. Both errors are computed for bothpossible mapping directions, i.e., A and B are interchanged in(60)-(63). To obtain the mean geometry for the one-way er-ror, the two geometries to be averaged need to have the sameglobal orientation. For this reason, their heart axes, as deter-mined in section 2.4.5, are aligned before averaging. After av-eraging the Euclidean coordinates, the mean geometries have to be remeshed, because the numerical computation of ventric-ular coordinates requires meshes of su ffi cient quality, which isnot guaranteed after moving the nodes. Linear interpolation isused to transfer the ground truth coordinates onto the remeshedgeometries. We use fTetWild (Hu et al., 2020) for remeshing.As an envelope size of 5 % of the mean edge length is used toapproximate the original mesh, the influence of remeshing onthe results can be considered negligible. For an example of theactual geometries and coordinates involved in the computationof the one-way error, the reader is referred to Fig. 24. The linearity of coordinates is important to preserve normal-ized distances. It is particularly relevant for transferring cardiacactivation times, where a shortening or lengthening in space canlead to artificial regions of slow or fast conduction, respectively. et al. (2021)
To evaluate the linearity of the rotational coordinate, we ex-tract contour lines at discrete isovalues of the apicobasal andtransmural coordinates and plot the normalized distance alongthese contour lines over the rotational coordinate. Ideally, theresult should be a diagonal line passing through (0 ,
0) and (1 , a , a (cid:48) ∈ (cid:110) , , . . . , (cid:111) (64) r , r (cid:48) ∈ (cid:110) , , . . . , (cid:111) (65) m , m (cid:48) ∈ (cid:110) , , . . . , (cid:111) (66)This yields 90 contour lines for the evaluation of the rotationaland 180 contour lines for the evaluation of the apicobasal linear-ity error. Linear interpolation was used to resample the normal-ized distance along the contour lines at 1000 equidistant valuesof the rotational (apicobasal) coordinate.
4. Results
Cobiveco was successfully and autonomously computed onall 1000 quasi-random instances of the SSM and all 36 patientgeometries, which demonstrates the robustness of the method-ology and its implementation.
Fig. 14 provides a visual comparison of Cobiveco and UVCfor all four coordinates on the mean shape of the SSM and twoexemplary patient geometries.As the mean shape has a very uniform wall thickness, the con-tour lines of the rotational and apicobasal coordinates appearequidistant for both methods, but artifacts at the discontinuitiesof the rotational coordinate can be seen for UVC (green circles).Patient 36 also has a relatively uniform wall thickness, but dif-ferences between both methods become more apparent. ForUVC, the distance between contour lines of the rotational co-ordinate increases near the septal junctions (magenta vs. cyancircle), which is not the case for Cobiveco.In patient 33, the di ff erences are most pronounced. While thecoordinates computed using Cobiveco still change very uni-formly in space, there are substantial distortions in the UVCcoordinates. The length of the segments between contour linesof the rotational coordinate changes up to four-fold betweenregions of small and large wall thickness. The apicobasal coor-dinate is also distributed very non-uniformly, indicating that thegeodesic approach to normalize the apicobasal Laplace solutiondoes not work reliably. In fact, a slight change of the geometrycan cause a di ff erent geodesic path between apex and base to become the shortest and therefore lead to an abrupt change ofthe apicobasal coordinate. Taking a closer look at the trans-mural coordinate within the LV shows that it changes muchfaster at the endocardium than it does at the epicardium becausethe width of the region between the two boundary surfaces in-creases with the circumference.If the coordinates always showed the same distortions for ev-ery geometry, this would only be a minor problem. However,comparing the rotational and apicobasal UVC coordinates forpatient 33 and the mean shape reveals that the same coordinatevalues can represent quite di ff erent anatomical regions (yellowstars). In contrast, the coordinates obtained using Cobiveco areconsistent across the geometries (green stars).For pictures showing Cobiveco and UVC coordinates on all 36patient geometries, the reader is referred to Fig. 22 and 23, re-spectively. The mapping errors as defined in section 3.3 were computedfor all patient geometries as well as both possible mapping di-rections. To condense the results, we averaged the error his-tograms across all geometries and both mapping directions.This leads to an equal weighting of errors for each case, in-dependent of the number of nodes in the respective mesh. Theaverage histograms are depicted in Fig. 15. Statistical measures(vertical lines) of the average histograms are given in Table 1.The two-way error shows a 3.5-fold improvement of the meanand the 99 th percentile is reduced even more. However, the me-dian is increased, which indicates that there are more small( < .
013 mm), but fewer large errors than for UVC. With amean value well below one mean edge length, our two-way er-rors for UVC are comparable to those in Bayer et al. (2018).The one-way error is more relevant in practice as it goes beyondevaluation of interpolation errors. Here, the error histogram de-cays much faster for Cobiveco and all statistical measures showa more than 4-fold improvement compared to UVC. In particu-lar, the mean one-way error is reduced from 7.1 to 1.5 mm andthe 99 th percentile is reduced from about 24 to 6 mm.Fig. 16 shows mapping errors for each individual patient ge-ometry. In all patients, the 99 th percentile of the two-way errorfor Cobiveco is below the mean edge length, which is not thecase for UVC (note the broken y-axis). For one-way errors, thelargest 99 th percentile in a single patient is about 8 mm for Co-biveco and 38 mm for UVC. To assess the spatial distributionof mapping errors, we visualized their mean across all patientson the mean shape of the SSM. To avoid artifacts due to spatialinterpolation, only errors directly available on the mean shapeof the SSM were taken into account for this purpose, i.e., onlyone mapping direction was included.The result in Fig. 17 (left) clearly shows that the two-way errorsfor UVC concentrate at discontinuities of the coordinates (com-pare with Fig. 1). Furthermore, there are large errors at the sin-gularities of the rotational coordinate. For Cobicevo, these er-rors are greatly reduced, because only the transventricular coor-dinate is discontinuous and the origin of the apicobasal coordi-nate coincides exactly with the rotational singularities. Choos-ing narrower colormap limits to visualize the two-way errors for . Schuler et al. (2021) 13 C ob i vec o U V C P a t i e n t C ob i vec o U V C M ea n s h a p e o f SS M C ob i vec o U V C P a t i e n t v , v ′ m , m ′ r , r ′ a , a ′ Fig. 14.
Visual comparison of coordinates ( v , m , r , a ) and ( v (cid:48) , m (cid:48) , r (cid:48) , a (cid:48) ) computed using Cobiveco and UVC, respectively, on the mean shape of the SSM and twoexemplary patient geometries. Green circles mark artifacts at discontinuities of r (cid:48) . Magenta and cyan circles mark regions of stretched and compressed coordinatevalues, respectively. Green and yellow stars mark an exemplary point with the coordinates (0 , , , ) for Cobiveco and UVC, respectively.4 S. Schuler et al. (2021) Table 1.
Summary of mapping errors (values of the vertical lines in Fig. 15).All values in mm.Error type Coordinate system Median Mean 90 th P. 99 th P.Two-way Cobiveco 0.013 0.038 0.084 0.385UVC 0.007 0.134 0.164 2.944Improvement factor 0.52
Cobiveco (Fig. 21) reveals the pattern of the isocurves used tocompute the apicobasal coordinate (Fig. 12, bottom-left). Thesemany non-zero, but still small errors explain the slight increasein the median two-way error observed for Cobiveco.For the one-way error (Fig. 17, right), the discontinuities andsingularities only play a minor role. It is dominated by incon-sistencies of the coordinates across di ff erent geometries, whichlead to inconsistent point correspondences. On average, thelargest one-way errors occur at the RV outflow tract for Co-biveco and at the apical region of the LV lateral wall for UVC.Nevertheless, absolute errors are much smaller for Cobiveco. The linearity of the rotational and apicobasal coordinate wasevaluated as described in section 3.4. As several thousand con-tour lines were extracted from the patient geometries, plottingthe normalized distance over the respective coordinate for eachindividual contour line would yield too many curves for visualinterpretation. Therefore, we created 2D histograms of thesecurves. The result can be seen in the first row of Fig. 18. Thesecond row shows the actual linearity error, i.e., the absolute de-viation from the black diagonal in the first row. Here, the meanand the standard deviation were computed across the di ff erentcontour lines, but not across the points along a contour line. Thethird row shows the dependence of the rotational (apicobasal)linearity error on the apicobasal (rotational) coordinate. Here,the mean and the standard deviation were computed along andacross all contour lines with the same apicobasal (rotational)isovalue. Table 2 summarizes the linearity errors using the max-ima of the mean curves in the second row of Fig. 18. In contrastto UVC, the apicobasal coordinate of Cobiveco shows almostperfect linearity. This is expected, as its computation is basedon isocontours of the other coordinates. For both methods, therotational linearity error is largest in the RV, but Cobiveco alsoshows a more than 4-fold improvement. As the circumferenceof the ventricles increases from apex to base, the (relative) ro-tational linearity error should decrease with the apicobasal co-ordinate. This can only be observed for Cobiveco. Table 2.
Summary of linearity errors (maximum of the curves in the secondrow of Fig. 18). All values in %.Coordinate Coordinate system RV mean (std) LV mean (std)Rotational Cobiveco 1.64 (1.59) 1.09 (1.25)UVC 7.11 (5.48) 5.04 (3.63)Improvement factor (3.44) (2.91)Apicobasal Cobiveco 0.65 (0.25) 0.51 (0.19)UVC 3.53 (2.99) 5.44 (4.05)Improvement factor (11.82) (21.65)
5. Application examples
One potential application of Cobiveco is the visualization ofcardiac data. Apart from the transfer of data from di ff erenthearts onto one common biventricular geometry for compara-tive visualization, the coordinates can also be used for a pro-jection of data onto a 2D representation of the 3D geometry.As a standardized way of visualization, we suggest to representthe surface of the biventricular myocardium using three polarprojections: One for the epicardium of both ventricles, one forthe RV endocardium, and one for the LV endocardium. Fig. 19shows an example for visualization of a geodesic distance fieldoriginating at the center of the RV septal surface. This examplewas chosen, as it allows a visual assessment of geometric dis-tortions caused by the projections. The main advantage of thepolar projections (lower half) is that the entire surface is visi-ble, whereas large regions remain obscured in the correspond-ing 3D views (upper half) even after individual rotation of thethree surfaces. Polar projections obtained using Cobiveco arean alternative to the method in Stoks et al. (2020), which usescylindrical coordinates to project a ventricular surface onto acone and then onto a circular disk and is limited to the epi-cardium or the LV endocardium only. Cobiveco also allows tocreate polar projections for any transmural layer between theendo- and epicardium. To obtain the projections, polar coordi-nates are computed for a cartesian grid with the desired targetresolution. The radial and angular polar coordinates are theninterpreted as apicobasal and rotational ventricular coordinates,respectively, and a mapping matrix is constructed as describedin section 2.5. We provide a function for computing polar pro-jections as part of the Cobiveco code. (cid:1)(cid:2)(cid:3)(cid:4)(cid:5)(cid:6)(cid:7)(cid:3)(cid:8)(cid:9) (cid:1)(cid:2) (cid:3)(cid:2) (cid:1)(cid:10)(cid:7)(cid:11)(cid:4)(cid:5)(cid:6)(cid:7)(cid:3)(cid:8)(cid:9) (cid:1)(cid:2) (cid:4)(cid:5)(cid:6)(cid:7)(cid:8)(cid:9) (cid:3)(cid:2) Geodesic distance (mm)
Fig. 19.
Visualization of a geodesic distance field using polar projections.
Top:
Original data on the epi- and endocardial surfaces.
Bottom : Corresponding po-lar plots with projected data.
Black dots : Apex at a = Black lines : Transven-tricular / septal junctions at r = r = / et al. (2021) 15 C ob i vec o U V C Twoway error (mm) P r obab ili t y Oneway error (mm) P r obab ili t y Twoway error (mm) P r obab ili t y Oneway error (mm) P r obab ili t y Fig. 15.
Average histograms of two-way errors ( left ) and one-way errors ( right ) evaluated for Cobiveco ( top ) and UVC ( bottom ). Histograms were averaged acrossall 36 patients and, for both types of errors, include both possible mapping directions. Each histogram contains about 50 M data points. C ob i vec o U V C Patient number T w o w a y e rr o r ( mm ) Patient number O ne w a y e rr o r ( mm )
25 1 6 11 16 21 26 31 36
Patient number T w o w a y e rr o r ( mm ) Patient number O ne w a y e rr o r ( mm ) Fig. 16.
Bar charts showing statistical measures of mapping errors for each individual patient. The data are the same as in Fig. 15. Patient 6 was excluded from theevaluation of UVC because the rotational and apicobasal UVC coordinates were too inconsistent to obtain a proper mean geometry for the one-way error.6 S. Schuler et al. (2021)
Mean of twoway errors (mm) C ob i vec o U V C Mean of oneway errors (mm)
Fig. 17.
Spatial distribution of the mean mapping errors across all patients visualized on the mean shape of the SSM. As a common geometry is needed to averageerrors across patients, only mapping sequences with respect to the mean shape of the SSM are included here. RV r N o r m a li z ed d i s t an c e LV r r L i nea r i t y e rr o r ( % ) r a L i nea r i t y e rr o r ( % ) a RV r' LV r' r' r' RV a N o r m a li z ed d i s t an c e LV a a L i nea r i t y e rr o r ( % ) a r L i nea r i t y e rr o r ( % ) r RV LV r' r'Sqrt of relative frequency of occurence Mean Mean + Standard Deviation
Cobiveco UVC Cobiveco UVCLinearity of rotational coordinate Linearity of apicobasal coordinate a' a'a' a'a' a'
Fig. 18.
Linearity of the rotational ( left ) and apicobasal ( right ) coordinate of Cobiveco and UVC. In the first row , the normalized distance along the contour linesis plotted over the coordinate to be assessed for linearity. A 2D histogram of curves corresponding to the individual contour lines is shown color-coded. Ideally,all points should lie on the black diagonal. The absolute deviation from the black diagonal (linearity error) is plotted in the second row . The third row shows thedependency of the linearity error on the respective other coordinate (apicobasal coordinate a or a (cid:48) for the rotational linearity error and rotational coordinate r or r (cid:48) for the apicobasal linearity error). r (cid:48) is only defined in the interval [0 , /
3] for UVC in the RV.. Schuler et al. (2021) 17
Another application is the integration of data from elec-troanatomical mapping and tomographic imaging. Fig. 20shows an example for the transfer of activation times recordedusing the CARTO mapping system (Biosense Webster, Inc.,Irvine, USA) onto the corresponding surfaces of a volume meshcreated from CT images. To compute Cobiveco, the endocar-dial surfaces from CARTO were first converted into a volumemesh (see Fig. 25 for a rule-based pipeline to create a vol-ume mesh from only endocardial surfaces). The coordinatesobtained on both geometries were then utilized to transfer theactivation times.In contrast to nearest-neighbor mapping (Duchateau et al.,2018; Graham et al., 2019) or other straightforward meth-ods (Cedilnik et al., 2018), Cobiveco allows a continuous andbijective mapping between geometries from both modalities.We believe that an unwanted smoothing of activation timesshould not motivate a discontinuous mapping between both ge-ometries (Duchateau et al., 2018) but should be addressed by anappropriate spatial upsampling on the source geometry.
CARTO mesh CT-derived mesh
Activation times (ms)r
Fig. 20.
Transfer of activation times recorded using CARTO.
Upper half : TheCARTO mesh is converted into a volume mesh and Cobiveco is computed forboth meshes.
Lower half : The coordinates are used to map the activation timesfrom the CARTO mesh to the endocardial surface of the CT-derived mesh.
6. Discussion
The evaluation of mapping and linearity errors showed thatCobiveco o ff ers a more consistent description of biventricularposition than UVC. These improvements are practically rele-vant. In the context of ECGI, for example, localization errorslie in the order of 10 to 30 mm (Potyagaylo et al., 2019; Gra-ham et al., 2020). The use of ventricular coordinates for vali-dation of non-invasive cardiac mapping or for machine learn-ing based approaches to this problem is only justified if er-rors due to the coordinates are substantially smaller. There-fore, especially the reduction of the one-way error from 7.1 to1.5 mm (mean) and from 24 to 6 mm (99 th percentile) is im-portant. A drawback of Cobiveco compared to UVC is the in-creased computational complexity. On a modern personal com-puter (8 × ffi ciency of our implementation could be im-proved through parallelized remeshing, parallelized isocontourextraction and more advanced preconditioning of linear systemsif computational e ff ort becomes crucial for certain use cases. Cobiveco has the following limitations: • The angle between local directions of the apicobasal andthe rotational coordinate can become small. This was es-pecially observed near the RV outflow tract (see patient 3in Fig. 22, for example) and might explain the larger one-way errors in this region (Fig. 17, top-right). As we de-cided for normalized and linear coordinates, the angle be-tween coordinate directions directly depends on the shapeof the geometry. • The transventricular coordinate remains discontinuous andis defined by a Laplace solution. Although the limitationsof Laplace solutions for this purpose are not as severe asfor the other coordinates, there might be more accurateways to separate both ventricles at the center of the sep-tal wall. • For all Laplace solutions, zero Neumann conditions areimposed at boundaries without Dirichlet conditions. Thisis reasonable for the ridge Laplace solution in (26) andfor the Laplace solutions used to obtain the tangent fieldsin (17) and (39). However, natural boundary conditionsas suggested in Stein et al. (2018) might be more appro-priate for the Laplacian extrapolation in (50) and for thetransventricular Laplace solution in (10). • The coordinate system does not cover myocardial bridgesbetween the atrioventricular valves and the outflow tracts.
7. Conclusion
We compared di ff erent approaches to define and computeventricular coordinates and developed Cobiveco, a consistentbiventricular coordinate system. Key novelties and improve-ments of Cobiveco are the symmetry of coordinate directions in et al. (2021) both ventricles and the definition of coordinates values based onthe normalized distance along bijective trajectories between twoboundaries, which can be computed by solving linear PDEs. Toavoid errors due to imprecise internal boundaries, we use im-plicit domain remeshing. The resulting coordinates are continu-ous (apart from the binary transventricular coordinate), normal-ized, and change linearly in space. To assess the consistencyof the coordinates across di ff erent geometries, a novel one-waymapping error was introduced. Evaluation on 36 patient ge-ometries showed a more than 4-fold reduction of mapping andlinearity errors compared to UVC. These improvements makeCobiveco an accurate analysis tool and a reliable building blockfor data-driven modeling of the cardiac ventricles. Acknowledgments
The authors would like to thank Olaf D¨ossel, AndreasWachter, Luca Azzolin and Gerald Moik for valuable discus-sions and feedback.This work received funding by the German Research Associ-ation (DFG) under grants DO 637 / / References
Bai, W., Shi, W., de Marvao, A., Dawes, T.J.W., O’Regan, D.P., Cook, S.A.,Rueckert, D., 2015. A bi-ventricular cardiac atlas built from 1000 + highresolution MR images of healthy subjects and an analysis of shape and mo-tion. Medical image analysis 26, 133–45. doi: .Bayer, J., Prassl, A.J., Pashaei, A., Gomez, J.F., Frontera, A., Neic, A., Plank,G., Vigmond, E.J., 2018. Universal ventricular coordinates: A genericframework for describing position within the heart and transferring data.Medical image analysis 45, 83–93. doi: .Cedilnik, N., Duchateau, J., Dubois, R., Sacher, F., Ja¨ıs, P., Cochet, H., Serme-sant, M., 2018. Fast personalized electrophysiological models from com-puted tomography images for ventricular tachycardia ablation planning. Eu-ropace : European pacing, arrhythmias, and cardiac electrophysiology :journal of the working groups on cardiac pacing, arrhythmias, and cardiaccellular electrophysiology of the European Society of Cardiology 20, iii94–iii101. doi: .Cerqueira, M.D., Weissman, N.J., Dilsizian, V., Jacobs, A.K., Kaul, S., Laskey,W.K., Pennell, D.J., Rumberger, J.A., Ryan, T., Verani, M.S., Schroeder,W., Schroeder, W., Martin, K., Lorensen, B., Schroeder, W., Schroeder, W.,Martin, K., Lorensen, B., 2002. Standardized myocardial segmentation andnomenclature for tomographic imaging of the heart. Circulation 105, 539–542. doi: .Chmelevsky, M., Budanova, M., Zubarev, S., Potyagaylo, D., Treshkur, T.,Lebedev, D., 2018. Clinical evaluation of noninvasive ECGI epi-endocardialmapping accuracy, in: Computing in Cardiology (CinC), pp. 1–4.Cluitmans, M., Brooks, D.H., MacLeod, R., D¨ossel, O., Guillem, M.S., vanDam, P.M., Svehlikova, J., He, B., Sapp, J., Wang, L., Bear, L., 2018.Validation and opportunities of electrocardiographic imaging: From tech-nical achievements to clinical applications. Frontiers in Physiology 9:1305.doi: .Costa, K.D., Hunter, P.J., Wayne, J.S., Waldman, L.K., Guccione, J.M., Mc-Culloch, A.D., 1996. A three-dimensional finite element method for largeelastic deformations of ventricular myocardium: II—prolate spheroidal co-ordinates. Journal of Biomechanical Engineering 118, 464–472. doi: .Crane, K., Weischedel, C., Wardetzky, M., 2013. Geodesics in heat. ACMTransactions on Graphics 32, 1–11. doi: . Dapogny, C., Dobrzynski, C., Frey, P., 2014. Three-dimensional adaptive do-main remeshing, implicit domain meshing, and applications to free and mov-ing boundary problems. Journal of Computational Physics 262, 358–378.doi: .Duchateau, J., Sacher, F., Pambrun, T., Derval, N., Chamorro-Servent, J., De-nis, A., Ploux, S., Hocini, M., Ja¨ıs, P., Bernus, O., Haissaguerre, M., Dubois,R., 2018. Performance and limitations of noninvasive cardiac activationmapping. Heart rhythm doi: .Geuzaine, C., Remacle, J.F., 2009. Gmsh: A 3-d finite element mesh generatorwith built-in pre- and post-processing facilities. International Journal for Nu-merical Methods in Engineering 79, 1309–1331. doi: .Graham, A.J., Orini, M., Zacur, E., Dhillon, G., Daw, H., Srinivasan, N.T.,Lane, J.D., Cambridge, A., Garcia, J., lly, N.J., Whittaker-Axon, S., Tag-gart, P., Lowe, M., Finlay, M., Earley, M.J., Chow, A., Sporton, S., Dhinoja,M., Schilling, R.J., Hunter, R.J., Lambiase, P.D., 2019. Simultaneous com-parison of electrocardiographic imaging and epicardial contact mapping instructural heart disease. Circulation. Arrhythmia and electrophysiology 12,e007120. doi: .Graham, A.J., Orini, M., Zacur, E., Dhillon, G., Daw, H., Srinivasan, N.T.,Martin, C., Lane, J., Mansell, J.S., Cambridge, A., Garcia, J., Pugliese, F.,Segal, O., Ahsan, S., Lowe, M., Finlay, M., Earley, M.J., Chow, A., Sporton,S., Dhinoja, M., Hunter, R.J., Schilling, R.J., Lambiase, P.D., 2020. Eval-uation of ECG imaging to map hemodynamically stable and unstable ven-tricular arrhythmias. Circulation: Arrhythmia and Electrophysiology 13.doi: .Hu, Y., Schneider, T., Wang, B., Zorin, D., Panozzo, D., 2020. Fast tetrahedralmeshing in the wild. ACM Trans. Graph. 39. URL: https://doi.org/10.1145/3386569.3392385 , doi: .Jacobson, A., 2013. Algorithms and Interfaces for Real-Time Deformationof 2D and 3D Shapes. Ph.D. thesis. ETH Zurich. Z¨urich. doi: .Jacobson, A., 2018. gptoolbox: Geometry processing toolbox.Jakob, W., Tarini, M., Panozzo, D., Sorkine-Hornung, O., 2015. Instant field-aligned meshes. ACM Transactions on Graphics 34, 1–15. doi: .de Marvao, A., Dawes, T.J.W., Shi, W., Minas, C., Keenan, N.G., Diamond,T., Durighel, G., Montana, G., Rueckert, D., Cook, S.A., O’Regan, D.P.,2014. Population-based studies of myocardial hypertrophy: high resolu-tion cardiovascular magnetic resonance atlases improve statistical power.Journal of cardiovascular magnetic resonance : o ffi cial journal of theSociety for Cardiovascular Magnetic Resonance 16, 16. doi: .Paun, B., Bijnens, B., Iles, T., Iaizzo, P.A., Butako ff , C., 2017. Patient inde-pendent representation of the detailed cardiac ventricular anatomy. Medicalimage analysis 35, 270–287. doi: .Potyagaylo, D., Chmelevsky, M., Budanova, M., Zubarev, S., Treshkur, T.,Lebedev, D., 2019. Combination of lead-field theory with cardiac vectordirection: ECG imaging of septal ventricular activation. Journal of Electro-cardiology 57, S40–S44. doi: .Reinsch, C.H., 1967. Smoothing by spline functions. Numerische Mathematik10, 177–183. doi: .Revishvili, A.S., Wissner, E., Lebedev, D.S., Lemes, C., Deiss, S., Met-zner, A., Kalinin, V.V., Sopov, O.V., Labartkava, E.Z., Kalinin, A.V.,Chmelevsky, M., Zubarev, S.V., Chaykovskaya, M.K., Tsiklauri, M.G.,Kuck, K.H., 2015. Validation of the mapping accuracy of a novel non-invasive epicardial and endocardial electrophysiology system. Europace :European Pacing, Arrhythmias, and Cardiac Electrophysiology : Journal ofthe Working Groups on Cardiac Pacing, Arrhythmias, and Cardiac CellularElectrophysiology of the European Society of Cardiology 17, 1282–1288.doi: .Roney, C.H., Pashaei, A., Meo, M., Dubois, R., Boyle, P.M., Trayanova, N.A.,Cochet, H., Niederer, S.A., Vigmond, E.J., 2019. Universal atrial coordi-nates applied to visualisation, registration and construction of patient spe-cific meshes. Medical image analysis 55, 65–75. doi: .Schroeder, W.J., Martin, K., Lorensen, W.E., Schroeder, W., Schroeder, W.,Martin, K., Lorensen, B., 2006. The Visualization Toolkit: An Object-Oriented Approach to 3D Graphics. Kitware Inc., Clifton Park, NY.Stein, O., Grinspun, E., Wardetzky, M., Jacobson, A., 2018. Natural boundaryconditions for smoothing in geometry processing. ACM Transactions onGraphics 37, 1–13. doi: .Stoks, J., Nguyen, U.C., Peeters, R., Volders, P.G.A., Cluitmans, M.J.M., 2020.. Schuler et al. (2021) 19An open-source algorithm for standardized bullseye visualization of high-resolution cardiac ventricular data: UNISYS, in: Computing in Cardiol-ogy. URL: .Yang, T., Yu, L., Jin, Q., Wu, L., He, B., 2018. Localization of origins ofpremature ventricular contraction by means of convolutional neural networkfrom 12-lead ECG. IEEE transactions on bio-medical engineering 65, 1662–1671. doi: .Yezzi, A.J., Prince, J.L., 2003. An eulerian PDE approach for computing tissuethickness. IEEE transactions on medical imaging 22, 1332–9. doi: .Zhou, S., Abdel Wahab, A., Sapp, J.L., Warren, J.W., Hor´aˇcek, B.M., 2019.Localization of ventricular activation origin from the 12-lead ECG: A com-parison of linear regression with non-linear methods of machine learn-ing. Annals of biomedical engineering 47, 403–412. doi: . Supplementary material
Table 3.
Tetrahedral mesh statistics of the mean shape of the SSM andthe 36 patient geometries. N : Number of nodes, l e : mean edge length, L : length along the long axis, W : width along the left-right axis. Minimum andmaximum values for the patient geometries are marked in bold.Mesh(Pat N (k) l e (mm) L (cm) W (cm) Mesh(Pat N (k) l e (mm) L (cm) W (cm)SSM 479 0.82 10.2 11.31 900 0.66 7.9 11.5 19 1103 0.91 9.9 13.52 825 0.69 7.5 11.2 20 1096 0.82 8.8 13.53 1029 0.71 8.6 11.1 21 1216 0.96 11.5 13.94 730 0.82 9.7 11.6 22 1146 0.96 11.7 14.55 1105 0.93 9.5 15.3 23 881 0.86 11.1 12.26
29 1327 0.76 9.4 11.412 872 1.01 11.8 14.7 30 934 1.01 12.1 14.813 912
31 1283 0.80 9.6 12.414 977 0.94 11.2 14.4 32 959 0.83 9.6 12.415 1105 0.83 9.9 12.7 33 972 0.81 10.4 11.616 862 0.88 10.3 12.8 34 1391 0.68 8.4 9.817 924 0.82 10.7 12.1 35 703 0.75 8.4 11.918 1268 0.81 10.3 12.2 36
Mean of two-way errors (mm)
Fig. 21.
Visualization of the mean two-way errors for Cobiveco as in Fig. 17,but with narrower colormap limits. The pattern of isocurves used to computethe apicobasal coordinate is visible.0 S. Schuler et al. (2021)
Fig.
22. Cobiveco computed for 36 patient geometries.
First rows : Rotational coordinate r . Second rows : apicobasal coordinate a .. Schuler et al. (2021) 21 Fig.
23. UVC computed for 36 patient geometries.
First rows : Rotational coordinate r (cid:48) . Second rows : apicobasal coordinate a (cid:48) .2 S. Schuler et al. (2021) Target (patient 33)Source (mean shape of SSM) Mean using Cobiveco Mean using UVC C ob i vec o U V C Ground truthTo be evaluated To be evaluated Ground truth r , r ′ a , a ′ Fig. 24.
Geometries and coordinates involved in the computation of the one-way error for patient 33.
Upper panel : Source, target and mean geometries.
Lowerpanel : Coordinates to be evaluated and ground truth coordinates on the mean geometry.
A B C D E
Distance to endocardium (mm) Laplace solution Distance to basal plane (mm) -20 200 0 1 -10 1000.5
Fig. 25.