Echocardiogram-based ventricular isogeometric cardiac analysis using multi-patch fitted NURBS

Monitoring the cardiac function often relies on non-invasive vital measurements, electrocardiograms, and low-resolution echocardiograms. The monitoring data can be augmented with numerical modeling results to support treatment-risk assessment. Often, medical images are not suitable for high-fidelity modeling due to their spatial sparsity and low resolution. In this work, we present a workflow that converts a limited number of two-dimensional (2D) echocardiogram images into analysis-suitable left-ventricle geometries using isogeometric analysis (IGA). The image data is fitted by reshaping a multi-patch NURBS template. This fitting procedure results in a smooth interpolation with a limited number of degrees of freedom. In addition, a coupled 3D-0D cardiac model is extended with a pre-stressing method to perform ventricular-cardiac-mechanic analyses based on in vivo images. The workflow is benchmarked using population-averaged imaging data, which demonstrates that both the shape of the left ventricle and its mechanical response are well predicted in sparse-data scenarios. The case of a clinical patient is considered to demonstrate the suitability of the workflow in a real-data setting.


Introduction
For decades, there has been an increasing trend in the use of patient-specific computational models in the clinical practice [1][2][3][4].Such models help to improve the understanding of biophysical phenomena and can provide decision support for clinicians [3,4].In the context of cardiac modeling, patient-specific computational models have been used to study a broad range of phenomena and diseases, such as valvular diseases [5,6], dyssynchronous heart failure [7][8][9], atrial arrhythmias [10][11][12], and ventricular arrhythmias under various conditions [2,[13][14][15][16].These models are most commonly based on high-quality image data (e.g., from magnetic resonance imaging (MRI) [7,10,15] or computer tomography (CT) [5,6]), but modeling based on sparse, irregular, or low-resolution image data (e.g., from echocardiograms [17]) is also possible.This work focuses on patient-specific analyses based on echocardiogram imaging data, which is typically only available in a limited number of cross-sectional planes (i.e., the data is sparse and irregular) in which only relatively large features can be identified properly (i.e., the resolution of the data is low).
The choice of a particular computational modeling technique is strongly influenced by the characteristics of the imaging data.Computational techniques to extract an analysis-suitable geometry from scan data can be classified as either direct segmentation techniques or as template-fitting techniques.In direct segmentation techniques, a mesh for the object of interest is obtained by converting imaging cells directly into computational elements.This typically requires scan data with sufficient resolution throughout the complete domain.The marching cubes method [18,19] -in which voxels are meshed based on an interpolation of the grayscale intensity -is the most prominent direct Preprint submitted to Computer Methods in Applied Mechanics and Engineering February 13, 2024 segmentation method in the context of the finite element method (FEM) [20][21][22].Scan-based immersed analysis methods [23] -in which the direct segmentation takes place on a background mesh that covers the scan domain -are a common state-of-the-art alternative.In template-fitting techniques, a reference geometry that captures the essential features of the object of interest is (gradually) deformed to match the scan data.This approach is less dependent on the resolution of the scan data, as the template can interpolate the geometry in regions where data is sparse.In the context of the FEM, typically the boundary of a template mesh (e.g., a surface triangulation (STL)) is deformed based on imaging data, after which the interior elements are updated in accordance with the boundary deformation to ensure sufficient mesh quality [24].Alternatively, template geometries based on splines can be used in the context of isogeometric analysis [25].While both fitting paradigms have their own merits, the suitability of template-fitting techniques for sparse and irregular data scenarios makes it the favorable approach in the context of this manuscript.In the anatomical modeling stage, a parametrized NURBS (Non-Uniform Rational B-Splines) template of the left ventricle, which is denoted by Ω lv templ , is deformed according to the sparse echocardiogram data, D. The resulting fitted geometry is deemed the image configuration, Ω lv img , that is subject to an image cavity pressure, p lv img .A rule-based fiber field, with direction e f img , is then constructed based on the image configuration.The geometrical properties, Ω lv templ , e f img , and p lv img , resulting from the anatomical modeling stage are then used as input for the cardiac modeling stage, which solves the relevant physics and consists of a three-dimensional (3D) mechanical component and a zero-dimensional (0D) circulatory component.An essential part of the cardiac modeling stage is the unloading of the image configuration -gradually decrease of p lv img to 0 -to retrieve the (stress-free) reference configuration, Ω lv 0 .

Anatomical model Cardiac model
The template-fitting technique proposed in this work builds on the mechanical cardiac modeling approach proposed by Willems et al. [26], which performs an isogeometric analysis (IGA) based on a multi-patch NURBS [27] spline model of the ventricles (illustrated in Fig. 1a).The main advantage of the proposed IGA approach -in which both the reference geometry and the deformations are interpolated using the same spline basis -is that the curved ventricles can be accurately represented using a minimal number of degrees of freedom (in the form of control points).Additional degrees of freedom are then only inserted in regions where this is needed to accurately mimic the mechanical response.The mechanical cardiac model of Ref. [26] has been benchmarked based on a state-of-the-art finite element solver and has been extended in Ref. [16] to incorporate the effect of infarct regions on the cardiac response, coupled with electrophysiology simulations.The IGA approach is demonstrated to make accurate predictions with a substantial reduction in the number of degrees of freedom compared to a traditional FE approach.Additionally, the absence of geometry clean-up and meshing steps in IGA makes it suitable for studying geometry variations.
Patient-specific analyses can be performed using the approach of Willems et al. [26], in the sense that the parameters of the spline model (e.g., cavity volume, wall thickness) can be retrieved from imaging data.Although the obtained ventricle models will follow the general characteristics of the scanned object, geometric details are lost in the parameter-retrieval step.In this regard, the proposed IGA approach does not fully exploit the geometric information encoded in the scan data.However, the spline-based ventricle model is an ideal starting point for a template-fitting technique.The main idea of this work is to extend the approach of Ref. [26] to take into account geometric details, as illustrated by the workflow in Fig. 1.While in the original approach the mechanical analysis is directly performed on the template, in this work we perform an intermediate step in which we deform the spline-based ventricle template based on the scan data using a template-fitting technique.This additional step will match the scanned object more precisely, generating more accurate patient-specific analysis results.
We herein assume that the echocardiogram data has been segmented, resulting in a three-dimensional point cloud (or polyline) representation of the boundary of the ventricle in each echocardiogram plane.In combination, these planar point clouds result in a sparse and irregular point cloud representation of the ventricle boundaries.A templatefitting technique based on a spline-based ventricle model is anticipated to be powerful in this setting, as the template can be fitted to the data points in regions where data is available (i.e., in the echocardiogram slices), whereas the splines will provide a smooth interpolation based on the template in the data scant regions (i.e., in between the slices).Although there exist many algorithms to fit splines to point clouds (see, e.g., Refs.[28][29][30] for overviews of various algorithms), these are often not directly applicable to the echocardiogram-based isogeometric analysis setting considered in this work.Specific challenges in this regard are the sparsity and irregularity of the point cloud and the use of a multi-patch NURBS object which is desired to retain smoothness across the patch interfaces.These challenges make techniques that approach the fitting as a point distance minimization problem promising in the context of the current work, since these allow for the simultaneous minimization of a distance measure and a smoothness measure [29].Our work incorporates various aspects established in the literature on such techniques, most notably the use of an iterative procedure for updating the points on the parametrized spline surface closest to the data points [31,32].Another aspect included in our work is the consideration of spline refinements to improve the quality of the fit [33].While the fundamental ingredients of the spline-fitting algorithm required in this work are established, an automatic algorithm capable of generating an analysis-suitable multi-patch NURBS object from the specific data set under consideration is not available.
A key contribution of the current work is the development of a fitting algorithm for multi-patch NURBS based on point distance minimization, suitable for sparse, irregular, and low-resolution echocardiogram imaging data.Novel features that improve the robustness and efficiency of the developed algorithm include the consideration of a convolution-based control point update scheme -which draws inspiration from scan-based immersed isogeometric analysis procedures [23] -and the use of a Voronoi tessellation to quantify the importance of data points relative to the fitting process.The developed fitting algorithm extends the patient-specific workflow of Willems et al. [26], allowing the IGA framework to exploit the full potential of the echocardiogram scan data.The extended workflow is tested using both dense and sparse data sets, considering both the geometry approximation and the mechanical cardiac response.With respect to the cardiac analysis, in this paper, we also extend the model of Ref. [26] to adjust for the fact that the imaging data is typically not obtained in an unloaded state.As part of this extension, the general pre-stressing algorithm proposed by Weisbecker et al. [34] is adapted to the considered isogeometric setting.We demonstrate the proposed patient-specific isogeometric analysis workflow in a real-data setting using a clinical patient.
This paper is outlined as follows.Section 2 commences with the introduction of the spline-based anatomical model.Section 3 then introduces the fitting algorithm, including examples to demonstrate its performance.Section 4 then reviews the cardiac model of Ref. [26] and introduces the extension to accommodate the pressure-loaded image configuration.The IGA workflow is then validated using dense point cloud data in Section 5, where the impact of data sparsity is investigated.An application scenario for a clinical patient is then studied in Section 6, after which conclusions and recommendations are presented in Section 7. described by a control net, consisting of i × j control points, P (P) i, j,2 , and corresponding weights, w (P) i, j,2 .For visibility, we only show the epicardial control net.

The spline-based anatomical model
The patient-specific isogeometric analysis workflow considered in this work requires the definition of a suitable spline-based computational domain.The multi-patch anatomical model of the idealized left ventricle developed in Willems et al. [26] serves as the template geometry in our patient-specific analyses.While the reader is referred to Ref. [26] for details, in this section we will review the nomenclature essential to this work.

The multi-patch NURBS template
We employ a common simplification of the left ventricle using a thick-walled ellipsoid that is truncated at the basal plane (Fig. 2 (a)).This idealized model assumes a smooth representation of the epi-and endocardium and neglects the papillary muscles and trabeculations inside the ventricles, which are responsible for opening the mitral and tricuspid valves.The idealized ventricle is represented by a multi-patch Non-Uniform Rational B-spline (NURBS) volume, which, in contrast to regular B-splines, enables the exact representation of conic sections.The physical geometry is obtained by mapping a multi-patch parameter domain to the physical domain, i.e., the truncated ellipsoid, using a geometrical map based on the NURBS basis functions, for which the control net provides the corresponding coefficients.This mapping is C 0 -continuous at the patch interfaces, in principle allowing for kinks in the geometry at the patch interfaces, but the control net is selected such that the geometry is (almost) smooth, i.e., G 1 -continuous.
We define the left ventricle template by the open bounded computational domain Ω lv ⊂ R 3 .The boundary of this domain comprises the epicardial (outer) surface, Γ epi , the endocardial (inner) surfaces, Γ endo , and the basal plane surface, Γ base .The closure of the computational domain is denoted by Ω lv , which is partition by n patch = 5 subdomains, Ω lv,(P) , referred to as patches, where P = {1, ..., n patch } is the patch index.The closed domain, Ω lv , is then defined as the union of the non-overlapping patches: The volumetric patches, Ω lv,(P) , are mapped from a parametric domain, Ω ⊂ R 3 , with coordinates (ξ, η, ζ), in accordance with where P (P) i, j,k ∈ R 3 are the patch control points.The NURBS basis functions are defined as where w (P) i, j,k are the control point weights and are tensor-product trivariate B-splines.The coordinate ζ corresponds to the through-the-thickness direction, with ζ = 0 corresponding to the inner surface, and ζ = 1 to the outer surface.The coordinates ξ and η correspond to the directions tangential these surfaces.The reader is referred to Ref. [28] for further details on the construction of B-spline bases.
In this work, for the left ventricle geometry parametrization second-order splines are used for the tangential directions and first-order basis functions are used for the through-the-thickness direction.In this setting, the complete volumetric control net is specified through the control nets of the inner and outer surfaces in accordance with where P (P) i, j,1 are the control point associated with the inner surface and P (P) i, j,2 those associated with the outer surface.From the perspective of the fitting algorithm this means that when the fitting is performed for the inner and outer surface individually, a complete parametrization of the ventricle is attained.In order to perform the isogeometric analysis, further refinement (knot insertion) and order elevation procedures may still be required.

The sparse data fitting algorithm
In this section, we present the developed multi-patch NURBS fitting algorithm suitable for sparse and irregular point clouds.Since the left ventricle geometry is completely defined through the inner and outer surfaces, we introduce the fitting procedure in the generic setting of fitting a multi-patch NURBS surface to a point cloud.To not overcomplicate notation, the surface to be fitted is denoted by Ω.In the context of the left ventricle, this means that Ω = Γ endo or Ω = Γ epi .The control net and NURBS basis functions corresponding to the surface Ω are here denoted by P (P)  i, j and R (P) i, j (ξ, η), respectively, simplifying the notation as used for the surface definitions (5).In Section 3.1 we introduce the point distance minimization problem underlying the fitting procedure, along with the developed iterative algorithm used to solve this problem.In Section 3.2 we then analyse the developed fitting algorithm using a well-defined test case.

The fitting algorithm
The fitting algorithm discussed in this section is developed to deform multi-patch surface topologies in 3D space, given a set of arbitrary (spatially sparse) data points.The data points can come from different sources, but in this manuscript we limit ourselves to 2D echocardiogram data (Figure 1).Such echocardiogram data relies on 2D slices or planes that, combined in 3D space, result in a sparse point cloud.We introduce the test case visualized in Figure 3 as an exemplary sparse point cloud fitting problem, and use it throughout this section to explain each algorithm component.The data points (indicated in blue) are defined in 3 different planes and are sampled from a cylindrical shape with unit radius.The goal of the fitting algorithm is to interpolate these data points using a multi-patch domain that is analysis-suitable.We do this by defining an initial template multi-patch geometry which is iteratively deformed such that it matches the data points and maintains smoothness, i.e., retains G 1 -continuity.
Let Ω be a closed multi-patch domain comprised of four patches as illustrated in Figure 3 (a).The domain consists of external boundaries {Γ v ∪ Γ h } ⊂ Ω and four internal boundaries, Γ int , i.e., interfaces.Given a set of n data arbitrary spatial data points, D = {d 1 , d 2 , ..., d n data }, where each item corresponds to the spatial coordinate of a data point, we require a solution to the minimization problem: Find the closed multi-patch domain Ω fit such that The functional F (•) to be minimized is the summation of the data fit error and the continuity/smoothness error across the interfaces.The constant λ introduces dimensional consistency and balances the importance of the error contributions [29].
The displacement error norm is an averaged discrete L 2 -norm that quantifies how close the data points are to the multi-patch surface and is defined as with the error vector, E i , visualized in Figure 3 (b), given by The projected data points, d proj i , are defined as the nearest-point-projection between the multi-patch surface and the data points (indicated in red in Figure 3 (b)), such that The continuity error norm, ε C , is quantified by the average angle between the normal vectors, n, across the interfaces, i.e., where the superscripts refer to the patch indices that share the considered interface, Γ int .Note that this error norm assumes that a maximum of two patches are allowed to share a single interface.Both error terms are chosen such that they represent a physically interpretable measure of the fitted surface, facilitating the user in choosing a suitable stopping criterion for the considered application.The algorithm aims to reduce both error contributions by changing the spatial coordinates of the control net, P i, j , following Equation (2), in an iterative sequential manner.This procedure, outlined in Algorithm 1, is graphically illustrated in Figure 4.The sequential algorithm comprises three main steps: (I) data point fitting of the individual patches (Figure 4 (a) and Lines 6-8 in Alg.1); (II) connecting the individual patches (Figure 4 (b) and Line 11); and (III) approximating G 1 inter-patch continuity (Figure 4 (c) and Line 13).Steps (I)-(III), which will be elaborated upon in the remainder of this section, are iterated, resulting in a gradual minimization of the functional in Equation (6).Each step is proceeded by a control net and surface update, (Lines 10, 12 and 14), where only the spatial coordinates of the control net, P i, j , are altered while the NURBS weights, w i, j , remain unchanged (visualized in Figure 4 (a)).Therefore, the set of control points that defines the multi-patch domain, P i, j = {P (1)  i, j , ..., P (P) i, j }, is updated every iteration (indexed by t) by a set of displacement vectors, u i, j = {u (1)  i, j , ..., u (P) i, j } for every patch P, as i, j = P (P) t i, j + C (u (P) i, j ), (11) where the P (P) t+1 i, j is the updated control net of patch P and C (•) is a user-defined constraint operator.This operator ensures that any computed displacement, u (P) i, j , adheres to user-defined constraints, e.g., is only allowed to move within a plane, along a line, or is completely fixed.Note that a relaxation factor, α relax , is used in Line 8 to control the computed displacement, u (P)  i, j , before updating the control net, cf.Equation (11).
Algorithm 1 Multi-patch template fitting algorithm.
Input:   The first step of the fitting algorithm involves displacing the multi-patch surface toward the data points.We do this for each patch separately, where we treat each patch as a separate surface that is not connected to adjacent patches, i.e., the patches are discontinuous.To approximate the control net displacement vectors u (P)  i, j of a single patch (Line 8 in Algorithm 1), we draw inspiration from the B-spline convolution operation considered in Ref. [35], which has favorable properties in terms of, e.g., the occurrence of oscillations.In the context of the setting considered here, the NURBS-based convolution operation is given by where the mismatch in data-point-to-surface is represented by the error vector, (Equation ( 8)).The error vector dictates the direction and distance the projected surface point, d The projected data points, d ), following Equation (2).The rational basis functions in Equation ( 12) are weighted by a mass m (P)  k .This mass represents the importance of the specific data point, d k , in relation to the fitting process and relative to other data points associated with the same patch.The data point mass is related to the surface area of Voronoi segments corresponding to the local projected data point, d local,(P) k , visualized in Figure 5 (a).These segments partition the parameter domain in regions that are closest to each of the local projected data points.In the proposed algorithm, the Voronoi segments are evaluated for each patch individually and are bounded by the local minima and maxima, i.e., ξ, η ∈ [0, 1] 2 .As a result, the Voronoi segments are not continuous across patch interfaces.Satisfying interface continuity of the Voronoi segments is conceptually possible, but constitutes a non-essential extension of the algorithm which is beyond the scope of this manuscript.Nonetheless, the projected data points that are located on a patch interface are still shared or seen by the patches connected to it.The Voronoi segments are defined in the parametric space, {ξ, η} ⊂ [0, 1] 2 , resulting in a total unit area A (P) = (P) k m (P) k = 1.Using the partition of unity property of NURBS, the area can also be expressed as a sum of control point masses, M (P) i, j , as provided that n (P) data 0 and with These masses indicate how important the control points are to the fitting process relative to one another and will be used in the continuity correction steps discussed below.In the case of an empty set of data points, n (P) data = 0, Equation ( 15) becomes a zero-valued matrix M (P) i, j = 0 i, j .A graphical representation of Equation ( 15) is visualized in Figure 5 (b), where the control point values are visualized as a heat map at the first iteration indicated for the third patch, Ω (P) .

3.1.2.
Step II: G 0 -patch continuity (Algorithm 1, Line 11) After displacing the individual patches following Equation ( 12), G 0 -patch continuity is not ensured and should be corrected for as is visualized in Figure 4 (a).Provided that we assume the patch interfaces to remain conforming, i.e., the knot vector at each interface side is identical, G 0 -continuity at the closed interface between two patches, Γ , is satisfied when for all control points on the interface, where K is a positive real number [36].It requires the spatial coordinates of the control points, P, at the shared interface to match identically, while the weights associated with these control points, w, may differ by a factor of K. Given that we keep the weights constant during the fitting process, we merely have to match the control points at the interfaces.However, one should consider the importance that the individual patches have on the data point fit, i.e., a simple average of the interface and vertex control points would hamper the fitting procedure if one of the patches does not contribute to the data fitting.This is the case for the first patch, Ω , visualized in Figure 3 (b), where no projected data points are attributed to this patch (indicated by the red dots).As a result, we perform a weighted average of the control points located at the internal interfaces and shared vertices cf.Equation (15), to obtain a G 0 -continuous multi-patch domain (visualized in Figure 4 (b)).
To perform the weighted average we first make a distinction between open interfaces and interface vertices, i.e., the corner points of patches on interfaces.Let χ int = {x k int | k ∈ I x int } be the unique set of interface vertices, where the set I x int comprises the set of patch indices adjacent to this vertex.For the example case in Figure 3 (a) we have I x int = {(1|2), (2|3), (3|4), (1|4), (1|2|3|4)}.Similarly, we define as the unique set of open multipatch interfaces, where the set I Γ int comprises the set of patch indices adjacent to this interface.For the example case in Figure 3 (a) we have I Γ int = {(1|2), (2|3), (3|4), (1|4)}.We define open interfaces such that they do not include their endpoints, which are interface vertices by definition, meaning Σ int ∩ χ int = ∅.By evaluating the displacement and distributed mass, Equation ( 12) and ( 15) respectively, for each patch that shares the control points located on the unique interfaces and vertices, Σ int , we compute the weighted average displacement such that, where the subscript λ indicates whether the quantity is evaluated at an open interface or interface vertex.This distinction between an open interface and an interface vertex is due to the difference in valency.The valency of an interface is typically limited to 2, with some exceptions, while the valency of interface vertices is usually higher [37].Because of this, open interfaces and interface vertices should be treated separately to ensure a correct weighted average is obtained following Equation (17).

3.1.3.
Step III: G 1 -patch continuity (Algorithm 1, Line 13) Following the correction on G 0 -patch continuity is the correction for G 1 -continuity, Figure 4 (c).This step aims to reduce the sharp G 0 edges at interfaces of interest and obtain a smooth G 1 representation of the geometry.These sharp edges are inherent to the multi-patch configuration due to the C 0 -continuity of the rational basis functions at the patch interfaces.Even if the parametrization is limited to C 0 -continuity at the patch edges, one can still obtain G 1continuity (or higher) as this is a geometric property rather than a parametrization argument.The geometric continuity, specifically G 1 , is quantified by the angle of the normal vectors (tangent planes) of each patch at the interface.While angle values should ideally be zero, angles with an absolute value below a specific threshold can be considered almost smooth (following Ref. [38] we use a threshold of 0.1 [deg]).We note that, in light of the developed fitting algorithm, the G 1 -continuity is not merely aesthetic, but plays an important role in maintaining the template shape during the fitting process, as will be discussed in the example below.
Satisfying G 1 -continuity is an active field of research in which the control points adjacent to and at the interface are positioned such that non-linear constraints ensuring the G 1 -continuity criterion are satisfied [36,[39][40][41].For arbitrary (complex) topologies, these non-linear constraints cannot always be satisfied, and as such often require degree elevation [42].In view of our analysis-suitable geometry, we do not desire excessive degree elevation or knot insertion since this would both defeat the purpose of fitting a coarse template on sparse data and increase computational load following the IGA paradigm.As a result, we propose a G 1 approximation procedure, which approximates G 1continuity at specific landmarks (set to be the Greville points).The number of landmarks is proportional to the number of control points, resulting in convergent G 1 -continuity behavior across the interface upon knot refinement, similar to Refs.[43,44].It should be mentioned that these geometry-based continuity methods are different from the continuity methods developed in the IGA framework: the Analysis-Suitable G 1 [45], the Approximate C 1 constructions [46], the D-patch [47], and the Almost-C 1 [48].Such continuity methods focus on the isogeometric function space rather than on the geometric property; see Ref. [49] for a comparison.Such alterations of the function space are not considered in this manuscript.
We consider again the multi-patch domain as shown in Figure 3, and provide a schematic overview of it in Figure 6 (a) with additional control points introduced for clarification.The control points that affect G 1 -continuity are located in the grey zone, which are either at the interface or adjacent to the interface (inner control points).The goal is to modify only the marked control points such that the multi-patch surface obtains globally approximate G 1continuity.The procedure of obtaining approximate G 1 -continuity across patch interfaces is explained using patches Ω (1) and Ω (2) as an example, Figure 6 (b).Unlike the open interfaces defined in the previous section, we consider only closed interfaces which are interfaces that include their endpoints (vertices), e.g., Γ int , visualized in Figure 3 (a).Let t (1) and t (2) be the tangent-boundary vectors and c (1) and c (2) be the cross-boundary vectors defined at the interface Γ (1|2) int of patches ( 1) and (2), respectively.Following the definition of G 0 -continuity in Equation ( 16), it holds that t (1|2) = t (1) = t (2) .G 1 -continuity is then guaranteed when the vectors, t (1|2) , c (1) , and c (2) , are linearly dependent across the entire closed interface, Γ (1|2) int .This means that the tangent space formed by, {c (1) , t (1|2) }, and, {c (2) , t (1|2) }, are both the same, i.e., the vectors are coplanar: det(c (1) , t (1|2) , c (2) ) = 0 [36].The procedure outlined in the following paragraph will therefore focus on modifying the specific control points such that this coplanarity condition is approximated.
Satisfying the coplanarity condition is achieved by first defining a common plane (C-plane) on the interface onto which the cross-and tangent-boundary vectors should be projected.In the general setting, this plane should be continuous along the interface, but satisfying continuity would put additional constraints on the spline degree [42].We therefore evaluate the G 1 -continuity at specific landmarks on the surface at which we minimize the G 1 -continuity jump in an iterative manner.We use the Greville abscissae [50] for the local coordinates of these landmarks because of their correlation to the curve.Let Ξi , i = 1, ..., n, and H j , j = 1, ..., m, be the Greville abscissae related to the knot vector Ξ and H.Following Equation (3), we define the global Greville points on the surface S (P) as In the remainder, G i, j , will be referred to as the Greville points, which are visualized in Figure 6 (b).For each row of Greville points that are located along the interface, denoted by the subscript i in Figure 6 (b), the same G 1 correction procedure is performed as explained in the next paragraph.Therefore, in the remainder, we only consider a single row and omit the row index i for notational brevity.Furthermore, Greville points located inside a patch, denoted by subscript j−1, are called inner Greville points.Greville points located on the interface Γ (1|2) int are denoted by, G (1|2) , and are related to the individual patch Greville points as, G (1|2) = G (1)  j = G (2) j , following the definition of G 0 -continuity in Equation (16).
The first step of the correction procedure is to construct two G (P) -planes at the inner Greville points of each patch, schematically visualized for the cross-section across the interface in Figure 7 (a) and graphically shown in (1) ∩ Ω (2) .The markers indicate the Greville points, G, associated with the G 1 -relevant control points in (a).For an interface with m rows of Greville points, indicated by the row index, we evaluate the cross-and tangent-boundary vectors, c and t respectively, at each Greville point G.These quantities are required to approximate G 1 -continuity across the interface. Ω Figure 7: Cross-sectional representation of two arbitrary NURBS patch surface, Ω (1) and Ω (2) , at an interface Greville point, G (1|2) , which can be, but is not a patch vertex in this case.(a) The first step involves the construction of two G (P) -planes at the inner Greville points, subscript j − 1, based on their normal vector, n (P) j−1 , and tangent-boundary vector, t (P) j−1 (out of plane component and not shown here).Two intersection points, A and C, are calculated based on the cross-boundary vectors at the interface, c (P) j .(b) The next step involves two line segments, BA and BC, on which two individual points, r (P) , are positioned based on the fitting relevance.The common C-plane is constructed based on Equation ( 23), of which the difference with the line segments' endpoints serves as the control point displacement correction, u (P) , cf.Equation (26). Figure 8 (b).The normal vector of these planes, n G (P) , is perpendicular to the normal, n (P) j−1 , and tangent-boundary vectors, t (P) j−1 , visualized in Figure 8 (a) and 8 (b).The G (P) -plane of patch P is then given as where Ω ( 2) Figure 8: Graphical overview of the G 1 -continuity correction steps visualized for patch Ω (1) and Ω (2) at Γ (1|2) int of the example case visualized in Figure 3. (a) The quantities of interest are first defined on the Greville points associated with the G 1 -relevant control points for a single interface-Greville point.(b) A combination of three planes is consecutively constructed, first the G (i) -planes, cf.Equation ( 19), and based on those the C-plane, cf.Equation ( 23).The C-plane serves as the common plane onto which the end-points of the line segments, BA and BC, are projected resulting in the control point displacement, u (P) , cf.Equation (26).Note that the result given in (c) is the final G 1 -correction result of the entire multi-patch surface, i.e., the visualized patches have also been influenced by the remaining patches, Ω (3) and Ω (4) , but are not visualized here.
Next, we define two line segments, BA and BC, such that with the lines, − − → BA and − − → BC, defined as The continuity plane, C-plane, is then defined as where The points r (1) and r (2) are located on the line segments, r (1) ⊂ BA and r (2) ⊂ BC, and their location on these segments depends on the data fitting information as defined in Equation ( 15), as follows r (1) = BA θ=θ (1) , with θ Considering the graphical representation in Figure 8, it is apparent that M (1) i, j = 0 i, j , since no data points are associated to patch Ω (1) .We then require that the control points that have the smallest influence on the data fit (low values of M i, j ) should be displaced more than control points that have a higher influence on the data fit.This ensures that the surface is minimally changed in areas that are important to the data fit during the continuity correction.We normalize the M (P) i, j weights for each line segment, BA and BC.We standardize the weights M (P) i, j individually for each line segment, namely BA and BC.This normalization is performed such that, in Equation ( 21), a value of θ = 0 signifies the highest significance for the interface (point B), while θ = 1 designates the inner control point (point A or C).For the special case that both segments yield θ = 0, resulting in r (1) = r (2) , we equate Equation (24) to the average of the normal vectors, 0.5(n (1)  j + n (2) j ).In the general setting, visualized in Figure 7 for two arbitrary NURBS surfaces, where both surfaces have equal importance for the data fitting, the points r (1) and r (1) are typically positioned close to the center of the line segments, Figure 7 (b).
The final step in satisfying G 1 -continuity is to displace the endpoints of the line segments BA and BC onto the constructed C-plane along a specified direction where subscript θ=0 and θ=1 are the parametric values corresponding to the endpoints of the line segments, cf.Equation (21).By applying the displacements of Equation ( 26) to the corresponding control points, P (1) j−1 , P (1|2) , and P (2) j−1 , the surfaces are displaced and G 1 -continuity is approximated.Equation ( 26 Here, the common C-plane is visualized onto which the endpoints of BA and BC are displaced.The resulting displacement of the control points yields an updated surface that approximates the G 1continuity across the patch interface, visualized in Figure 8 (c).However, the resulting multi-patch domain does not guarantee a continuous G 1 -continuity across all interfaces after a single G 1 correction update.Several iterations may be required before reaching a converged G 1 approximation due to the spline bases' (degree, knot refinement) inability to satisfy exact continuous G 1 -continuity at the interfaces.Nonetheless, the continuity is further improved upon knot insertion or spline degree elevation (the latter will not be considered in this manuscript).The proposed correction scheme can therefore be considered as an iterative and convergent scheme to achieve G 1 -continuity.Furthermore, the correction scheme is not unconditionally stable, often requiring constraints (Equation ( 11)) and relaxation on the computed displacement cf.Equation (26).A qualitative analysis of the continuity correction in combination with the data fit is examined in the next section.

Analysis of the fitting algorithm
The fitting algorithm outlined in Algorithm 1 is applied to the example case visualized in Figure 3.We apply constraints to the computed displacement during each fitting iteration to ensure stable convergence toward a physically valid final shape.Constraints ensure that the final multi-patch fit is regular, i.e., no singular control points [37], and that the continuity correction behaves stable for any arbitrary problem setting.To ensure this, we constrain all patch edges, i.e., the external boundaries and interfaces of the multi-patch domain in Figure 3, using the constraint operator, C (•) cf.Equation (11), such that These constraints ensure that the fitted geometry maintains its initial rectangular shape as good as possible.The inner control points of each patch are not constrained, which, for this test case, does not impede the stability of the fitting algorithm.In the analyses of the left ventricle in Section 5 we adopt a slightly different set of constraints due to the problem setting discussed there, which emphasizes the problem-specific nature of the constraints.The algorithm is run for t max = 400 fit iterations, during which a total of 3 uniform knot refinements are allowed once the gradient of the displacement error is below a specified tolerance, The point-tosurface projection, Line 3 in Algorithm 1, is calculated each fit iteration while the continuity-correction step is iterated  In Figure 9 (b) we plot the normalized displacement error with respect to the first iteration and the absolute average continuity error.The trace of the displacement error shows an exponential decay before reaching a plateau for each uniform refinement level.The plateau is caused by the inability of the initial non-rational B-spline multi-patch surface to precisely describe the cylindrical shape that the data points are based on [28].Only after uniformly inserting additional knots (Line 17), the initial surface can achieve < 0.1% displacement error.The continuity error of the multi-patch surface shows a different behavior, varying between refinement levels.We attribute this difference to the addition of knots: once the surface is uniformly refined it can locally adjust the shape of the geometry according to the data points.This local change in geometry introduces geometric gradients that impede the continuity correction as it becomes more difficult to satisfy G 1 -continuity globally.Additional knot refinements can relax the conditions needed to satisfy the global continuity conditions using our continuity correction scheme, see Section 3. The results presented in Figure 9 are based on an initial template with quadratic B-splines and no internal knots, i.e., a single element per patch, as visualized in Figure 3 (a).The initial knot vector of the multi-patch surface influences the final shape of the fitted geometry.To illustrate this, we use the knot vector from the fourth (final) surface in Figure 9 (a) as the starting knot vector for the initial template before running the fitting algorithm.The results are visualized in Figure 10, where a distinct difference is observed.Based on Equation ( 12) only the control points that have support on the projected data points are displaced, and the remaining control points are unchanged.As a result, only a small portion of the multi-patch domain is displaced, while the rest adheres to the initial template shape.This comparison shows the importance of the initial number of control points of the chosen template.A large number of control points leads to a minimal change in the geometry in areas where there is no available data to be fitted, while a small number of control points ensures a higher degree of change.In the case of the left ventricle, we minimize the initial number of control points for the template.This can only be done if the minimal number of control points can describe the desired template shape, i.e., an ellipsoid in this case, with sufficient accuracy.This emphasizes the strong benefit of splines to describe a geometry with a relatively low number of degrees of freedom (control points) over conventional polygonal techniques.

The isogeometric cardiac model
The echocardiogram-based electromechanical cardiac model considered in this work builds on the isogeometric analysis (IGA) approach proposed by Willems et al. [26].For self-containedness, in Section 4.1 we first briefly review the model of Ref. [26].In Section 4.2 we then discuss the model extensions required to accommodate the pressure-loaded image configuration.

The isogeometric cardiac model
In Section 4.1.1we first introduce the governing equations for the cardiac model, after which the solution techniques -in particular the spatial discretization using isogeometric analysis -are discussed in Section 4.1.2.We refer the reader to Willems et al. [26] for a detailed exposition of the model and the required solution techniques.

Governing equations
Our cardiac model consists of two main components, viz. a three-dimensional continuum model mimicking the electromechanical behavior of the left ventricle and a zero-dimensional lumped parameter model describing the circulatory system.These two model components are coupled through the left ventricle cavity volume and pressure.Below we discuss the most important aspects of the cardiac model, a comprehensive overview of which is presented in Box 1.
To accommodate the large deformations and deformation gradients that are physically essential, a total Lagrangian formulation is used for the three-dimensional electromechanical model.To define the total Lagrange formulation, a distinction is made between the undeformed configuration, Ω lv 0 , and the current configuration, Ω lv (Figure 11).The position of a material point in the undeformed configuration is denoted by x 0 i and that of a material point in the deformed configuration by x i .The total displacement of a material point is denoted by u i = x i − x 0 i .In the discussion in this section, we assume the undeformed configuration to be known, which allows us to present the cardiac model without the need to refer to the intermediate image configuration, Ω lv img (Figure 11).The incorporation of the image configuration will be discussed in Section 4.2.
The mechanical equilibrium equations are presented in Box1.1a-Box1.1g.Ignoring inertia effects and body forces, the mechanical equilibrium is described in the current configuration, Ω lv , by means of the Cauchy stress tensor σ i j (Box1.1a).On the endocardium, Γ endo , a pressure load p lv is applied (Box1.1b).On the basal plane, Γ base , the deformation of the ventricle is constrained in the normal direction (Box1.1c).The two remaining rigid body translations are constrained by prescribing the average of the in-plane components of the displacement field in the basal plane to be zero (Box1.1d-Box1.1e)and the rigid body rotation is constrained by prescribing the average curl of the displacement field to be zero (Box1.1f).These rigid body constraints are defined in the current configuration (see Appendix A).For the initial condition of the mechanical equilibrium equations a zero displacement field is considered (Box1.1g).
In the total Lagrangian formulation it is most convenient to express the constitutive behavior in terms of the Green-Lagrange strain tensor E i j (Box1.3)-which is directly related to the deformation gradient F i j (Box1.2) -and the second Piola-Kirchhoff stress tensor S i j .The Cauchy stress σ i j is related to the second Piola-Kirchhoff stress by Box1.5a.The second Piola-Kirchhoff stress in our model consists of two components, a passive component, S pas i j , corresponding to the hyperelastic behavior of the myocardium, and an active component, S act i j , corresponding to the force generated by the sarcomeres.
The passive component of the stress is expressed in terms of the invariants of the Green-Lagrange strain tensor I 1 , I 2 and I 3 (Box1.4a-Box1.4c)and the quasi-invariant I 4 (Box1.4d)representing the Green-Lagrange fiber strain.The fiber direction in the image configuration is determined using the model proposed in Ref. [51] and mapped toward the unloaded configuration as e f0 i .The passive stress component is described by the Fung-type constitutive model proposed by Bovendeerd et al. [52] (Box1.5b).This model assumes a transverse isotropic material that is parametrized by the elastic modulus of the myocardium, C, and the bulk modulus, κ.The quadratic strain form Q (Box1.5d)depends on the strain invariants I 1 , I 2 and I 4 , with a 1 , a 2 and a 3 being material constants.
The active component of the stress is described through the phenomenological model of Refs.[52,53], which assumes the contraction of the myocytes to be initiated simultaneously at time T act .In this model, a sarcomere of length l s is regarded as a serial connection of a contractile element with length l c , and an elastic element with length l s − l c .The stretch of the sarcomere is related to the Green-Lagrange strain in the fiber direction (Box1.5e).The length of the contractile element is governed by a contractile dynamics relation (Box1.6a),where v 0 represents the unloaded shortening velocity.Before activation, the contractile length is equal to the length of the sarcomere (Box1.6b).The stiffness of the elastic element of the sarcomere model is denoted by E a , and the stress generated by this element acts in the direction of the fibers in the undeformed configuration (Box1.5c).The rise and decay of the myocyte contraction are governed by the twitch function, f twitch , which depends on the time since activation t a = t − T act .The influence of the contractile length on the active stress component is prescribed by the isometric contraction function, f iso .For details regarding f twitch and f iso , including the used parameters, the reader is referred to Ref. [26].
The zero-dimensional model of the circulatory system is illustrated in Figure 12.In essence, the circulatory model provides a (differential) relation between the pressure in the left ventricle, p lv , the volume of the left ventricle cavity, V lv , and the arterial and venous pressure p = {p art , p ven }.The circulatory model is parametrized by the compliances, C art and C ven , the reference volumes, V art ref0 and V ven ref0 , and the resistances R art , R ven and R per .The volume of the left ventricle depends on the solution to the mechanical model (Box1.7b)and the mechanical model, in turn, depends on the pressure in the ventricle through the boundary condition (Box1.1b).This makes the system coupled in two directions.

Solution techniques
The formulation summarized in Box 1 is discretized in space using isogeometric analysis.In this discretization technique, the displacement field and contractile sarcomere length field are approximated by where N u i j are vector-valued B-spline basis functions for the displacement field and N l c i are scalar-valued B-spline basis functions for the contractile sarcomere length.The coefficients corresponding to the basis functions are denoted by ũi for the displacement field, and by lc i for the contractile sarcomere length.In this work, we consider maximum regularity cubic (p = 3) B-splines for the displacement the contractile length field.This choice of B-spline basis functions makes the discretized displacement and sarcomere field C 2 -continuous in the interior of the patches and C 0 -continuous across patch boundaries.
We consider a Bubnov-Galerkin discretization of the weak form of Box 1, which is expressed completely on the reference configuration.The pressure boundary condition is weakly enforced (through the weak form), while the normal displacement boundary condition on the basal plane is enforced strongly (by constraining the relevant control point displacements).The rigid body constraints are enforced through scalar Lagrange multipliers.The weak forms and their matrix representations resulting from this discretization are elaborated in Willems et al. [26].
Figure 12: Schematic overview of the 0D circulatory system coupled to the 3D LV domain, Ω lv .The closed-loop lumped parameter model represents a simplification of the human body's circulatory system.The circulatory system is modeled as two Windkessel compartments in series, characterized by a compliance, C i , and a volume at zero reference pressure, V i ref0 , with i ∈ {art, ven} (art=arterial and ven=venous).The compartments and ventricles are linked via resistances, R i (green), through which arterial (i = art), venous (i = ven), and peripheral (i = per) volume flow occurs, q i (orange), caused by the compartment pressure, p i (blue).The two heart valves: the mitral valve (MV) and aortic valve (AV), are idealized to open and close instantaneously and prevent back-flow.The ventricle and compartment pressures, illustrated in blue, are the quantities that are solved in the circulatory model.
Box 1 Overview of the electromechanical cardiac model based on Ref. [26].

Equilibrium equation:
) Deformation gradient: Green-Lagrange strain tensor: Invariants of the strain tensor: in Ω lv (Box1.4b) in Ω lv (Box1.4c) Constitutive relation: in Ω lv 0 (Box1.5d) Contractile dynamics: Circulatory model: For the temporal discretization of the time-dependent terms emanating from the contractile dynamics model and the circulatory system model we consider a backward Euler scheme.This results in a nonlinear system of equations for each time step, which is solved using Newton-Raphson iterations.To optimize computational effort while retaining stability, an adaptive time stepper is used.This time stepper scales the time step based on the required number of Newton iterations.

The pressure-loaded image configuration
A fundamental aspect in scan-based cardiac modeling is that the images are generally not obtained in an unloaded state.To take this effect into account, we extend the modeling approach of Ref. [26] so that it can accommodate a loaded scan configuration.The fundamental idea behind the proposed extension is that all computations are performed on the isogeometric anatomical model obtained by the fitting procedure introduced in Section 3, i.e., the image configuration.The rationale behind this is that the quality of the isogeometric geometry resulting from the fitting procedure can be assured, contrasting the case in which computations are performed on a geometric model deflated toward the unloaded configuration.In Section 4.2.1 we detail how the loaded image configuration is accommodated in the cardiac model.In Section 4.2.2 an isogeometric analysis version of the general pre-stress algorithm of Weisbecker et al. [34] is proposed to evaluate the mapping between the unloaded configuration and the image configuration.

The decomposition of displacements and displacement gradients
In our model, we make the distinction between three configurations, as schematically shown in Figure 11 for the ventricle.We denote the position of the material point in the unloaded configuration Ω lv 0 by x 0 i .When loaded by a cavity pressure p lv img , the ventricle is deformed toward the image configuration, Ω lv img .The displacement corresponding to this image configuration is denoted by u img i = x img i − x 0 i .When loaded with an arbitrary pressure, p lv = p lv img + p lv arb , the material is deformed toward the current configuration, in which the total displacement of the material points is given by where u arb i is the arbitrary displacement relative to the image configuration as a result of the arbitrary loading p lv arb .In the remainder, we will omit the superscript 'lv', e.g., p = p lv , for clarity but also to emphasize the generality of the GPA.
In our formulation, we assume that the image displacement, u img i , has been determined (in our case using the algorithm discussed in Section 4.2.2).The displacement field corresponding to an arbitrary pressure loading is then determined by the model elaborated in Section 4.1, but with the total displacement field in accordance with equation (30).To evaluate the model it is necessary to determine the total deformation gradient.From the above definition of the configurations it follows that x i = x 0 i +u img i +u arb i , from which the total deformation gradient can be determined as where The decomposition (31) conveys that when the image deformation, u img i (or the image deformation gradient, F img i j ) is known, the total deformation gradient F i j can be evaluated.Using the definitions in Box 1, the constitutive behavior can then be evaluated for an arbitrary loading condition.
In our formulation, we perform all operations on the image configuration.This implies that the gradient of u arb i with respect to the image configuration can directly be evaluated and that the corresponding deformation gradient F arb i j can be used to evaluate gradients and integrals in the deformed configuration in accordance with where f is an arbitrary function.
The image displacement, u img i , is also defined on the image configuration, such that the image deformation gradient in (32) cannot be evaluated directly.To evaluate this expression, the chain rule is applied to yield from which it follows that This reformulated expression for the (inverse of) the image deformation gradient can be evaluated since the gradient of the image displacement is now defined with respect to the image configuration.

The isogeometric general pre-stressing algorithm
To evaluate the cardiac model outlined in Section 4.1 in the case of a loaded image configuration, the image displacement u img (or the image deformation gradient field F img ) must be known.To determine this image displacement, we adapt the general pre-stressing algorithm (GPA) of Weisbecker et al. [34] to our isogeometric analysis setting.
The concept behind the GPA is to find the image deformation gradient field, F img , for which the passive response in the image configuration, x img , is in equilibrium with the applied pressure, p img .To formalize this concept, we express the passive model, i.e., the mechanical model comprised by the equilibrium equations (Box1.1a-Box1.1g)and the passive component of the constitutive relation (Box1.5b), in abstract form as This abstract form expresses that, given an image deformation gradient field F img and a pressure load p, the total deformation gradient field F can be computed.In the case that the pressure is equal to that of the image state, i.e., p = p img , the total deformation gradient field should be equal to the deformation gradient field of the image itself, that is This equilibrium condition for the image configuration is equivalent to stating that the arbitrary deformation gradient, F arb , in the decomposition ( 31) is equal to the identity tensor, and, following the definition (32) that the arbitrary displacement field vanishes, i.e., u arb = u − u img = 0. Conceptually, the nonlinear equation ( 37) can be solved for the image deformation gradient field F img , given a certain p img .The original algorithm proposed in Ref. [34] solves equation (37) by gradually increasing the image pressure up to p img , using (37) as an update rule for the image deformation gradient field, i.e., F k+1 img = M passive (F k img , p k ), where k is the iteration number.When the image pressure is reached, i.e., p k = p img , equation ( 37) is used to perform load-increment-free fixed point iterations, i.e., F k+1 img = M passive (F k img , p img ).For many (non-linear) problems, including the cases considered in this work, the gradual increase of the pressure before performing the fixed point iterations is required to ensure convergence.
To employ the general pre-stressing algorithm in the context of the considered isogeometric cardiac model, the original algorithm of Ref. [34] has been adapted.The proposed isogeometric analysis GPA is outlined in Algorithm 2. Compared to the original algorithm, the following modifications have been made: • In principle, in the GPA the image deformation field u img does not need to be determined explicitly, as the model can be evaluated directly based on the image deformation gradient field, F img .We find it practical, however, to explicitly compute the image deformation field, u img , on the same spline basis as used for the discretization of the displacement field, i.e., equation ( 29), as Algorithm 2 Isogeometric general pre-stressing algorithm Input: p img , n fixed # Image pressure, max free-load iterations Initialize: u 0 img = 0, p 0 = 0, ∆p 0 = ∆p max , k = 0 1: while (k − n load ) ≤ n fixed do 2: end if 18: end while where ũimg,i are the control point image displacements.This discrete displacement image field is then computed by the projection of equation (35) in the Frobenius norm as In this projection (Line 4), the constraints are applied in the same way as considered for the model discretization in Section 4.1.
Our choice to explicitly compute the image displacement field on the spline basis is practically motivated.In isogeometric analysis, the control point displacements are the natural way of storing the displacement field, in the sense that also the image geometry is represented by a control net.This means that an isogeometric analysis code will have the ability to interpolate the image displacement field in the same way as that it interpolates the geometry.The spline-based interpolation allows for the computation of the image deformation gradient field in the interior integration points (as used for the evaluation of the internal forces), on the endocardium (as used for the external forces) and on the basal plane (where Lagrange multiplier constraints are applied).It would be possible to store the image deformation gradient field in different quadrature point sets for the various contributions to the system of equations, but such a representation would not be standard.The representation in terms of the control point displacements standardizes the data format, which is beneficial from a computational workflow point of view.Although the spline interpolation (38) introduces an additional approximation step compared to the direct usage of the image deformation gradient in the quadrature points, this approximation places the image displacement in the same discrete function space as the arbitrary displacement field u arb .This means that the image displacement field is approximated in the same way, and with the same approximation properties, as the increments with respect to the image configuration.
• Due to the large deformations and high degree of non-linearity in the constitutive behavior of the considered isogeometric cardiac model, gradually increasing the pressure toward the image pressure is a necessity.To optimize the number of load steps required to reach the image pressure, an adaptive pressure step size is employed.
The pressure increment is adjusted according to the number of Newton iterations, n Newton , required to solve the isogeometric passive model up to a specified tolerance (Lines 6-7), where n Newton target is the desired number of Newton iterations.This adaptive load step procedure is initialized by prescribing ∆p 0 = ∆p max , and a load increment is redone with an adapted step size if the Newton solver does not converge within n Newton max iterations.
Convergence of the algorithm is measured via the increments in the image displacement (Line 9) as where use is made of the L 2 -norm over the image configuration.The image displacement is assumed to be converged when the full image pressure load is applied and the increment measure ( 40) is below a specified tolerance.The number of steps required to reach the pressure load is denoted by n load , and the number of load-increment-free iterations by n fixed .The influence of these numerical parameters is studied for benchmark problems in Appendix B.

Benchmarking the patient-specific workflow
To benchmark the developed patient-specific isogeometric analysis workflow (Figure 1), in this section we consider its application to a standard population-averaged data set (Cardiac Atlas Project [54]).This data set, based on cardiac MRI, provides a dense point cloud representation of the left ventricle at end-diastole (end filling and prior to contraction).We use the dense point cloud data as a reference case, and use it to construct point cloud representations of six synthetic echocardiogram slices: the apical 2-3-and 4-chamber views (AP2CH, AP3CH, and AP4CH) and the apical-, mid-and basal-parasternal-short-axis views (A-PSAX, M-PSAX, and B-PSAX) [55].The choice of synthetic slices is limited to and based on the available patient data that will be discussed in Chapter 6.However, additional views can be considered, based on a specific clinical echocardiogram protocol of interest.The various point sets are shown in Figure 13.The workflow is applied to the dense reference data (9744 points, labeled "dense"), to the sparse synthetic echocardiogram data using all six slices (575 points, labeled "spr.6"), and to the sparse synthetic echocardiogram data using only the AP2CH, AP4CH, and M-PSAX slices (344 points, labeled "spr.3"), which are the slices considered in the next chapter.We aim to quantify the effects of data sparsity by comparing the echocardiogram fitting results (Section 5.1) and the cardiac mechanics results (Section 5.2) to the results obtained using the dense reference case.

Fitting results
The ellipsoidal multi-patch NURBS template discussed in Section 2 is used as the initial condition of the fitting algorithm for all cases.The position and dimensions of the template are based on the stress-free geometry as discussed in Willems et al. [26] and based on Bovendeerd et al. [56].The template is visualized in Figure 13 at both the epicardium and the endocardium surfaces.The developed fitting algorithm is applied to the epicardium and endocardium separately, after which a volumetric multi-patch NURBS parametrization is obtained by through-the-thickness interpolation; see Section 2.1.The control points at the base, ∂Γ base , are only allowed to displace in a predefined plane.This plane corresponds to the XY-plane containing the average spatial position of the data points that are located at the base.Data points positioned above this plane are neglected as the fitting error cannot be minimized for these points due to the constraint.We use the same settings for the fitting algorithm as for the example case considered in Section 3.2, with the exception that t max = 500 iterations are considered, the allowed gradient of the displacement error is set to ∆ε D = 5 • 10 −3 , and the last 20 iterations are only used to correct for continuity.We note that the settings of the algorithm have not been optimized for the case under consideration and that the algorithm is robust in the sense that the results do not change erratically when varying the settings.Optimization of the settings, which is beyond the scope of the current work, can yield further improvements, but will not fundamentally alter the results.
In Figure 14 we trace the normalized displacement error (left axis, black lines) and continuity error ε C (right axis, blue lines) for the endocardium (Figure 16 (a)) and the epicardium (Figure 16 (b)).We adopt the displacement and continuity errors following Equations ( 7) and ( 10), where we normalize the displacement error by the initial displacement error at the first iteration, i.e., ε D /ε t=0 D .For all three cases, the error behaves similarly as for the example case considered in Section 3.2.The normalized displacement error decreases gradually, reaching a plateau when the best possible fit to the data is obtained.Upon refinement of the multi-patch NURBS (solid circle markers), the error decreases further.Note that because the refinement criterion is dependent on the gradient of the displacement error, refinements occur at different moments for the different cases.After 500 iterations and two refinement steps, the normalized displacement error has reduced to less than 3.5% percent of the initial template fit error for all cases.The remaining error is a result of the refined template not being able to capture the strongly curved geometry at the apex.For all cases, the continuity error (i.e., the average angle between interfaces) is within 0.55 [deg].On account of the larger curvature near the apex of the endocardium compared to the epicardium, the continuity error is observed to be larger there.
When comparing the errors for the different cases, the normalized displacement error is observed to decrease with the number of data points.The reason is two-fold: First, a reduction in data points reduces the fitting constraints, which ultimately improves the domain interpolation of the data points.Second, the initial displacement error, ε t=0 D , of the template is relatively higher compared to cases with a greater number of data points.This difference arises from the initial template and point cloud distribution.Provided that the initial template is already close to the data points (in regions away from the apex), a relatively lower initial error is computed for cases with an increased number of data points, which results from the averaging of the displacement error with respect to the number of data points, cf.Equation (7).Nonetheless, all three cases show a converged normalized displacement error of which the relative difference is negligible.Unlike the normalized displacement error, the continuity error is not directly influenced by the point cloud distribution.Consequently, there is a minimal difference in its convergence behavior across the various cases.
Figure 15 compares the final volumetric results -where the endocardial and epicardial boundaries are combined as explained in Section 2.1 -of the fitting algorithm for the three point-set cases.Visually, the results for the various cases are indiscernible.When considering the cavity and wall volumes, it is observed that for the six-slice case (Figure 15 (b)) these volumes deviate from the dense reference case by 0.2%.For the three-slice case (Figure 15 (c)) Circulatory R art 1.00 × 10 −2 kPa s ml −1 system R per 1.20 × 10 −1 kPa s ml −1 R ven 2.00 × 10 −3 kPa s ml −1 C art 2.50 × 10 1 ml kPa −1 C ven 6.00 × 10 2 ml kPa −1 V art ref0 5.30 × 10 2 ml V ven ref0 3.10 × 10 3 ml V tot 5.00 × 10 3 ml the uncertainties related to the scanning and segmentation process.In this light, for the considered setting we do not expect the errors related to the fitting procedure to be detrimental to the analysis workflow.

Cardiac mechanics results
Using the fitted image geometries obtained for the three cases we now evaluate the cardiac response using the model of Section 4.1.For all cases, we adopt parameter settings nearly identical to those employed in the model validation study by Willems et al. [26], with the only deviations being in the reference active stress, T 0 , and the reference volumes of the circulatory model, V art ref0 , V ven ref0 .These parameter values were slightly altered such that a physiological hemodynamic response was obtained.An overview is given in Table 1.The reader is referred to Ref. [26] for details regarding the choice of parameters.The final refinement of the fitting procedure, with two uniform refinement steps relative to the initial template, is directly used for the cardiac analysis.The IGA cardiac model is constructed using the Nutils Python toolbox [57], where we use maximum regularity cubic B-splines for the displacement and contractile length field, resulting in 3885 and 1351 degrees of freedom respectively.Based on the mesh convergence study presented in Ref. [26] this discretization setting is expected to yield accurate results.This expectation will be confirmed by consideration of a further uniform refinement at the end of this section.
In the first stage of the cardiac analysis, the unloaded configuration is determined using the algorithm of Section 4.2.The end-diastolic image pressure is assumed as 1600 [Pa] or 12 [mmHg].An initial pressure increment of ∆p 0 = 35 [Pa] is used for the isogeometric GPA algorithm.The pressure step is scaled using n Newton target = 7 and the iterations are terminated after n fixed = 20 load-increment-free iterations.The convergence plots for the GPA are shown in Figure 17.For all cases, the adaptive pressure load stepper -characterized by the nonlinear slope of the normalized pressure curve in Figure 17 (a) -reaches the image pressure after 27 iterations, marked by the solid spheres.The cavity-to-wall volume ratio, shown in Figure 17 (a), decreases from approximately 95% in the image configuration, to approximately 40% in the unloaded configuration.After the pressure is fixed at the image pressure, this ratio is virtually constant.From the relative displacement increment shown in Figure 17 (b), it is observed that the error decreases linearly while the pressure is fixed at the image pressure, showing the same behavior as for the GPA benchmark problems considered in Appendix B. The deflation procedure is further illustrated in Figure 18, which shows the deflation process for the dense point cloud case.Although the deflated geometry as shown in this figure is not directly used in the analysis, i.e., all operations are performed on the image geometry resulting from the fitting procedure, this plot clearly illustrates how the ventricle is deflated.It is particularly noteworthy to observe the torsional rotation that takes place during deflation on account of the fiber distribution.
With the image deformation state determined through the isogeometric GPA, the time-dependent cardiac model can be evaluated.We use an initial time step size of ∆t = 2 [ms], which is scaled by a factor of 1  2 in the case of a The cardiac-cycle duration is fixed at 800 [ms] of which five consecutive cycles are computed to achieve a cyclic hemodynamic steady state.The hemodynamic response of the fifth cycle is shown in Figure 19, showing very similar results for all cases.The minor differences between the cases can be explained directly from the differences in the fitting results.The maximum cavity volume loop is directly related to the cavity volume of the image state, which, as indicated in Figure 15, is smallest for the six-slice case, largest for the three-slice case, and with the dense point cloud   19) and (b) end-systole (circle marker in Figure 19).reference case in between.The deformation field obtained for the dense point cloud case is illustrated in Figure 20 at end-diastole and at end-systole (corresponding to the markers in Figure 19).The torsional motion of the ventricle is clearly observable from the mesh lines.The displacement results for the six-and three-slice cases are visually indiscernible from the dense case and have therefore not been included here.
To assess the influence of the mesh size, in Figure 21 the results obtained on the mesh considered above, with 5241 degrees of freedom in total, is compared to results with an additional uniform mesh refinement, which results in 22489 degrees of freedom.The obtained pressure-volume loops for the two considered meshes show very close resemblance.When considering local quantities, such as the myofiber strain in specific material points, minor differences are observed.These observations regarding the mesh convergence behavior are in agreement with the detailed mesh study reported in Ref. [26], to which the reader is referred for details.
As for the fitting results discussed above, the global cardiac response results convey that the data sparsity only marginally affects the results when compared to the dense point cloud reference case.With the settings of this benchmark study being chosen to be representative of real-world cases, including the way in which the synthetic echocardiogram data is being generated, the obtained results give confidence that the proposed workflow can yield reliable results in sparse data scenarios.We note that the observations presented here based on the global hemodynamic response do not necessarily carry over to other, possibly more local, quantities of interest.

Application to echocardiagram data
The patient-specific isogeometric analysis workflow -which has been benchmarked based on the Cardiac Atlas Project [54] -is designed to be applicable to real patient data.To demonstrate this applicability, in this section, we perform a cardiac analysis on the echocardiogram data of a clinical patient.
The scan data is obtained using a Philips iE33 ultrasound machine, with a spatial resolution of 640 × 480 and a temporal resolution of 37 frames.Three echocardiogram slices have been obtained, viz.apical 2-and 4-chamber views (AP2CH, AP4CH), and the parasternal-short-axis view (PSAX).The scanning results for each slice are stored in a DICOM [58] file, from which the image showing the end-diastolic configuration is selected based on the closure of the mitral valve.The endo-and epicardium are then segmented manually, by selecting points on these surfaces, which was validated by a cardiologist (Figure 22 (a)).The echocardiogram slices are combined in 3D space by positioning the AP2CH and AP4CH center at the base and with an angle difference of 60 [deg] [55].The PSAX view is then fitted based on these apical views in height and rotation, resulting in a three-dimensional point cloud, as shown in Figure 22 (b).
To initialize the fitting procedure, a truncated ellipsoid template with identical dimensions as mentioned in Section 5.1 was scaled by a factor of 0.8 and positioned below the mitral valve (Figure 22 (c)).The fitting algorithm is then executed with the same parameter settings as considered above but with t max = 220 fit iterations and 20 G 1continuity correction iterations for each fit iteration.These changes in the settings are necessary to ensure that the fitting algorithm did not over-correct for G 1 -continuity, due to the strong curvature in the provided data point set.The resulting fit -in which the point distance error is ≈ 7% compared to the initial template and the continuity error is 0.7 [deg] for the endocardium and 0.3 [deg] for the epicardium -is shown in Figure 22 (d).
An isogeometric cardiac analysis is then performed directly on the mesh obtained from the fitting procedure, using Nutils [57], which was refined 2 times compared to the initial template.This multi-patch NURBS geometry results in a total of 3885 degrees of freedom for the (cubic) displacement field and 1351 degrees of freedom for the (cubic) contractile length field.For the analysis, the same parameters are used as for these considered in Section 5, with the exception that the general pre-stressing algorithm iterations are terminated after 6 load-increment-free iterations to determine the pre-stresses in the image configuration.
The displacement field obtained from the cardiac analysis is illustrated in the deformed configuration at endsystole in Figure 22 (e).As for the benchmark cases considered above, the torsional deformation mode is clearly observed.Compared to the Cardiac Atlas displacement results, visualized in Figure 20 (b), this patient-specific case exhibits a similar deformation field and deformation magnitude.From the isogeometric analysis we also obtain the hemodynamic response, which is shown in Figure 22 (f), with the end diastole and end systole instance marked by a square and circle respectively.The deformation field and pressure-volume relation are here presented as typical results generated by the isogeometric analysis, but evidently, the full time-dependent mechanical and hemodynamic response is available for analysis.

Concluding remarks
In this work, we have presented a patient-specific cardiac modeling workflow based on the isogeometric analysis (IGA) paradigm.The proposed workflow has two distinctive features, viz.: (i) Patient-specific geometries can be constructed from sparse and irregular data, for example obtained from echocardiogram slices; (ii) The geometry construction operations and the physical analysis are seamlessly integrated, not requiring the use of intermediate geometry clean-up and meshing operations.
The proposed workflow builds on the IGA-based cardiac-mechanics model of Willems et al. [26], from which it also derives its distinctive features.When fitted to sparse data points, the multi-patch NURBS representation of the left ventricle is able to smoothly interpolate the geometry in the regions where data is scant.This same multi-patch NURBS can then be used directly as a basis for a Galerkin method, in this case discretizing the mechanical cardiac response of the left ventricle.To realize this functionality, in this work, we have presented two methodological innovations, viz.: (i) A multi-patch NURBS fitting algorithm for sparse data; (ii) An extension of the computational model of Ref. [26] enabling it to accommodate in vivo images.
The developed fitting algorithm is capable of dealing with sparse data by assigning data point weights to the various control points.In doing so, the data is closely fitted by the splines, while the splines interpolate the geometry in regions where data is absent.The considered multi-patch NURBS template in principle allows for the occurrence of kinks across patch interfaces, which often (and certainly also in this work) is undesirable.The fitting algorithm minimizes the occurrence of such kinks by amending the global error representing the distance of the geometry to the data points with a term that penalizes non-smooth patch interfaces.In the algorithm, these two error terms are then minimized sequentially through control net manipulations.
Objects observed in in vivo images are typically not in an unloaded state.To accommodate the loaded image state in our analysis workflow, the model of Ref. [26] has been amended with a multiplicative decomposition of the deformation gradient.In this decomposition, the image state is encoded in the deformation gradient field related to the mapping between the (unknown) unloaded configuration and the image configuration.The advantage of this approach is that all computational operations are performed on the image geometry, for which the mesh quality can be assured through the developed fitting algorithm.To determine the component of the deformation gradient associated with the image loading, the general pre-stressing algorithm originally proposed by Weisbecker et al. [34] has been adapted to the isogeometric analysis setting.The most prominent alteration to the algorithm pertains to the representation of the image deformation gradient field by an auxiliary displacement field.This representation by means of control point displacements makes the algorithm compatible with data standards in IGA, thereby promoting integration in existing implementations.
A benchmark study based on left-ventricle data from the Cardiac Atlas Project [54] has demonstrated the capabilities of the proposed workflow.In this study, synthetic sparse echocardiogram data is extracted from a dense reference point cloud, allowing us to assess the effects of data sparsity.The fitting algorithm is demonstrated to maintain a high degree of accuracy when considering a realistic selection of echocardiogram slices.When considering the hemodynamic response computed through the IGA model, only minor deviations from the reference solution are observed.These benchmarking results provide confidence in the suitability of the workflow for real data scenarios.To demonstrate the operation of the workflow in such a scenario, the case of a clinical patient is considered.While a detailed analysis of this specific scenario is beyond the scope of the current work, it demonstrates the robustness of the workflow.
Our patient-specific isogeometric analysis workflow has been developed in the context of cardiac modeling but is to a large extent generic.Specifically, the fitting algorithm is suitable to generate analysis-suitable splines for isogeometric analysis based on general point cloud data, both dense and sparse.The developed workflow is deterministic in nature, i.e., all quantities are assumed to be known precisely.The reality is, however, that uncertainties are present in all steps of the workflow.Depending on the case under consideration, the quantification of such uncertainties may be of interest [59][60][61].In order to open the doors to uncertainty quantification, we envision the need for the development of a reduced-order modeling approach, which would, e.g., allow for the application of stochastic sampling procedures.

Appendix A. Constraining of the rigid body modes
In our cardiac model, the normal component of the displacement vector is constrained at the basal plane (Box1.1c).Since the normal to the basal plane coincides with the x 3 direction (with x = (x, y, z)) this constraint can be enforced strongly by setting the u 3 displacement field component of the basal plane control points to zero.This constraint removes the rigid body motion in the x 3 direction and the rigid body rotations around the x 1 and x 2 axes.A common strategy to remove the remaining rigid body modes is to strongly constrain the in-plane displacements in one point (in the basal plane), making this a point around which the body can rotate.The rotational rigid body mode is then removed by considering another point in the basal plane and constraining the displacement component perpendicular to the vector between this point and the rotation point.In the patient-specific isogeometric analysis setting considered in this work, this method of constraining is not practical.Directly constraining the displacement in a point is most practical in the patch vertices, where the B-splines are interpolatory.To constrain the displacement in an arbitrary point linear constraints can be used, but, depending on the considered software framework, such constraints are not always practical from an implementation point of view.When using the patch vertices as constraining points, these are in general not positioned in such a way that the rotational rigid body mode can be constrained by a single component of the displacement field.Instead of resorting to linear constraints to solve this problem, in this work, we opt to constrain the rigid body modes by enforcing the conditions Box1.1d, Box1.1e and Box1.1f using Lagrange multipliers.
While the average displacement constraints (Box1.1dand Box1.1e) are standard, the constraint on the average curl of the displacement field (Box1.1f) is not.To explain how this constraint removes the rigid body rotation, we consider a deformed equilibrium state in the current configuration, x, and apply an arbitrary rigid body transformation to obtain a new configuration (see Figure A.23) where b 1 and b 2 are arbitrary translations in the x 1 and x 2 direction, and θ is an arbitrary angle of rotation around the x 3 axis.In the absence of the rigid body constraints (Box1.1d-Box1.1f),the transformed configuration, x ′ , satisfies the equilibrium conditions for any choice of the translations and rotation, making the problem ill-posed.
To demonstrate how the rigid body constraints (Box1.1d-Box1.1f)result in a unique solution, equation (A.1) is expressed in terms of the displacement field as  In our analysis, we only consider the bottom part of the domain, prescribing symmetry conditions (i.e., zero normal displacement and zero non-normal traction) at the mid-plane.The domain is discretized using 5 patches with 2 × 2 × 2 elements each.Cubic NURBS are used to parametrize the geometry, resulting in a control net with 135 control points.The boundary conditions are identical to Box1.1c, Box1.1d,Box1.1e, and Box1.1f, which are applied to the Γ mid .In contrast to the nearly incompressible rod benchmark, this simulation setup does not require the use of a line-search algorithm.It is observed that for all refinement levels, this error decreases under the application of the load and during the first load-free iterations.After approximately 17 steps (7 load-free steps) the errors reach a stable plateau.This plateau corresponds to the discretization error, as the displacement field required to reach the undeformed state cannot be represented exactly by the NURBS basis.Evidently, as the mesh is refined, the exact solution can be represented more accurately.As a consequence, the value at which the error reaches the plateau decreases under mesh refinement.The relative image displacement increment is shown in Figure B. 28 (b).As for the rod benchmark, the two stages of the GPA algorithm can be observed.In contrast to the volume ratio error, the increments continue to decrease once the discretization accuracy is reached.The reason for this is that the increments are measured relative to the discrete solution, meaning that it will vanish once the discrete solution is attained.As for the rod benchmark, the relative displacement increment provides an adequate metric to assess the convergence of the GPA.

Figure 1 :
Figure 1: The IGA workflow proposed in this work consists of two main stages: the anatomical modeling stage and the cardiac modeling stage.In the anatomical modeling stage, a parametrized NURBS (Non-Uniform Rational B-Splines) template of the left ventricle, which is denoted by Ω lv templ , is deformed according to the sparse echocardiogram data, D. The resulting fitted geometry is deemed the image configuration, Ω lv img , that is subject to an image cavity pressure, p lv img .A rule-based fiber field, with direction e f img , is then constructed based on the image configuration.The geometrical properties, Ω lvtempl , e f img , and p lv img , resulting from the anatomical modeling stage are then used as input for the cardiac modeling stage, which solves the relevant physics and consists of a three-dimensional (3D) mechanical component and a zero-dimensional (0D) circulatory component.An essential part of the cardiac modeling stage is the unloading of the image configuration -gradually decrease of p lv img to 0 -to retrieve the (stress-free) reference configuration, Ω lv 0 .

Figure 3 :
Figure 3: (a) Graphical representation of the considered test-case with Ω the closed initial multi-patch domain, i.e., the template, and d i ∈ D the data points to which the template is to be fitted, marked in blue.The multi-patch domain consists of four closed patches, Ω (P) ⊂ Ω for P = {1, 2, 3, 4}, connected via their interfaces, Γ (i| j) int = Ω (i) ∩ Ω ( j) with i j.The external closed boundaries are denoted by, {Γ v ∪ Γ h } ⊂ Ω.The control net vertices located on the interfaces are denoted as, x int , and represent the endpoints of a closed interface, Γ int .(b) Visualization of the nearest-point-projected data points, d proj i , on the multi-patch domain, indicated in red.The difference is given by the error vector, E i .

Figure 4 :
Figure 4: Graphical overview of the three main steps of the multi-patch surface fitting algorithm: (a) Fitting the individual patches to the data points; (b) Correcting for G 0 -patch-interface continuity, i.e., connecting the patches; (c) Correcting for G 1 -patch-interface continuity, i.e., smoothing the multi-patch geometry at the patch interfaces.

⊂ R 2 ,
proj,(P) k , should be displaced.The NURBS basis functions R (P) i, j are evaluated in the local data points, d local,(P) k .These local data points are defined in the (ξ, η)-parameter domain, d local,(P) k and are obtained by solving d local,(P) k = arg min dlocal ∈ [0,1] 2 proj,(P) k ⊂ R 3 , are related to the local data points by d proj,(P) k = S (P) (d local,(P) k

Figure 6 :
Figure 6: (a) Schematic overview of the multi-patch domain, Ω, in which 16 control points (cps) per patch are visualized for illustration purposes only.Control points located within the grey area are considered G 1 -relevant cps and comprise vertices, interface cps, and inner cps (indicated in teal).(b) Schematic illustration of the relevant quantities of interest at the interface, Γ(1|2) int = Ω (1) ∩ Ω(2) .The markers indicate the Greville points, ) is best explained schematically by Figure 7 (b) and graphically by Figure 8 (b).

Figure 9 :
Figure 9: (a) Overview of the different multi-patch configurations at different fitting iterations with different knot vectors or refinement levels.The color bar represents the angular jump, continuity error ε C , at the interfaces.(b) Trace of the displacement error ε D and continuity error ε C for a fixed number of iterations.The vertical dashed lines indicate the instance at which a uniform knot refinement is applied to the multi-patch surface, with the encircled numbers, 1 at t = 80, 2 at t = 140, 3 at t = 192, 4 at t = 400, corresponding to the configurations in (a).Uniform knot insertion ensures both a decrease in displacement and continuity error.
1.3.Both error quantities show a converged solution upon reaching 400 iterations, at which the displacement error has a value of 1.91 • 10 −3 and the continuity error 0.02 [deg].The latter is expected to be slightly above 0 [deg] since the presented correction scheme in Section 3.1.3requires further refinements to make the error vanish.

Figure 10 :
Figure 10: Comparison of two fitted multi-patch surfaces after 400 fit iterations with different initial knot vectors.(a) Multi-patch surface with initial 0 × 0 inner knots per patch and a final 7 × 7 inner knots per patch after uniform knot insertion during the fitting process.(b) Multi-patch surface with an initial and final 7 × 7 inner knots per patch, without any uniform knot insertion during the fitting process.

Figure 11 :
Figure 11: Schematic overview of the computational domain in (a) the undeformed configuration, Ω lv 0 ,(b) the image configuration, Ω lv img , and (c) the current configuration Ω lv , corresponding to an arbitrary deformation.

Figure 13 :
Figure 13: Dense point cloud data of a left ventricle at end-diastole obtained from the Cardiac Atlas Project [54] at (a) the epicardium and (b) the endocardium with the initial template visualized in the center.(c) Point cloud representations of synthetic echocardiograms generated from the dense point cloud in six slices.The longitudinal positions of the PSAX views are approximately based on 25%, 50%, and 75% of the total length of the ventricle.The ellipsoidal template used as the initial condition for the fitting algorithm is shown at (a) the epicardium and (b) the endocardium.

Figure 14 :
Figure 14: Convergence of the displacement error (left axis, black lines) and continuity error (right axis, blue lines) for (a) the endocardium and (b) the epicardium.The results obtained using different point sets based on the Cardiac Atlas Project are compared.

Figure 17 :Figure 18 :
Figure 17: Convergence plots of the isogeometric general pre-stressing algorithm: (a) the cavity-to-wall volume ratio, and (b) the relative displacement increments.The results obtained using different point sets based on the Cardiac Atlas Project are compared.

Figure 19 :
Figure 19: Hemodynamic response of the left ventricle with the parameters as in Table 1: (a) left-ventricle, arterial and venous pressure signals, and (b) left-ventricle cavity pressure-volume loop.The results obtained using different point sets based on the Cardiac Atlas Project are compared.

Figure 21 :
Figure 21: Comparison of (a) the pressure-volume loop and (b) the myocardium strain computed in four points at the equatorial plane (z = 0 in Figure 13 (c)), using a discretization Ω lv dense (ref = 0) with 5241 degrees of freedom and a uniform refinement, Ω lv dense (ref = 1), with 22489 degrees of freedom.

Figure 22 :
Figure 22: Overview of the cardiac IGA workflow applied to real patient data.(a) Three views of the 2D echocardiogram data, which is segmented for both the endo-(orange) and epicardium (blue) at end-diastole.(b) Combined echocardiogram slices in 3D, resulting in a sparse point cloud (fitted left ventricle already visualized).(c) Sparse point cloud resulting from the slice combination with the initial left ventricle template positioned in the center and slightly below the mitral valve (base).(d) Result of the fitting algorithm after 3 uniform refinements showing indiscernible visual agreement with the data points.(e) Displacement magnitude field resulting from the IGA cardiac simulation at end-systole (noticing the torsional movement).(f) The hemodynamic response corresponding to the cardiac simulation, indicating the end-diastole (square marker) and end-systole (circular marker) cardiac instances.
Figure A.23: Schematic illustration of an arbitrary domain that undergoes pure rotation around its global origin with an angle θ.(a) The resulting configuration experiences a displacement, u ′ , consisting of a translation, R(θ)x c − x c , of the center, x c , and a rotation, R(θ)x 0 − x 0 around this center, cf.Equation (A.1).(b) Removing the translational component, cf.Equation (A.3), ensures that the local rotation around the center is left, which should be removed cf.Equation (A.5).

Figure B. 25 :
Figure B.25: GPA iteration results for the axially loaded rod benchmark: (a) applied image force; (b) height of the unloaded geometry.The results for different load application steps are compared.

Figure B. 26 :
Figure B.26: Convergence metrics of the GPA iterations for the axially loaded rod benchmark: (a) absolute error in the height of the unloaded geometry; (b) relative increment of the displacement field.The results for different load application steps are compared.

Figure B. 27
shows the GPA results for various uniform refinements of the mesh, where the case of n ref = 1 corresponds to the 135 control points with which the geometry is parametrized.For all refinement cases, the image pressure is applied in 10 steps (Figure B.27 (a)).When considering the ratio between the cavity volume and the wall volume in Figure B.27 (b), it is observed that it converges to a value that is visually indiscernible from the exact value of 1 7 .The error in the volume ratio between the cavity and the wall is inspected closer in Figure B.28 (a).

Figure B. 28 :
Figure B.28: Convergence metrics of the GPA iterations for the pressure-loaded thick-walled sphere benchmark: (a) absolute error in the cavityto-wall volume ratio of the unloaded geometry; (b) relative increment of the displacement field.The results for different refinement levels are considered.
∆p max ← min ∆p max , p img − p k# Limit the maximum allowable step size by the image pressure

Table 1 :
[26]meters of the constitutive model components and the circulatory system model, taken from Willems et al.[26].