Comparison of interpolation methods of predominant cardiomyocyte orientation from in vivo and ex vivo cardiac diffusion tensor imaging data.

Cardiac electrophysiology and cardiac mechanics both depend on the average cardiomyocyte long-axis orientation. In the realm of personalized medicine, knowledge of the patient-specific changes in cardiac microstructure plays a crucial role. Patient-specific computational modelling has emerged as a tool to better understand disease progression. In vivo cardiac diffusion tensor imaging (cDTI) is a vital tool to non-destructively measure the average cardiomyocyte long-axis orientation in the heart. However, cDTI suffers from long scan times, rendering volumetric, high-resolution acquisitions challenging. Consequently, interpolation techniques are needed to populate bio-mechanical models with patient-specific average cardiomyocyte long-axis orientations. In this work, we compare five interpolation techniques applied to in vivo and ex vivo porcine input data. We compare two tensor interpolation approaches, one rule-based approximation, and two data-driven, low-rank models. We demonstrate the advantage of tensor interpolation techniques, resulting in lower interpolation errors than do low-rank models and rule-based methods adapted to cDTI data. In an ex vivo comparison, we study the influence of three imaging parameters that can be traded off against acquisition time: in-plane resolution, signal to noise ratio, and number of acquired short-axis imaging slices.

average cardiomyocyte long‐axis orientation

cardiac diffusion tensor imaging

cone of uncertainty

fractional anisotropy

heat flux coordinate

tensor interpolation using heat flux coordinates

tensor interpolation using prolate spheroidal coordinates

left ventricular

proper generalized decomposition

proper orthogonal decomposition

prolate spheroidal coordinate

signal to noise ratio

singular value decomposition

INTRODUCTION

Cardiovascular diseases are the leading cause of death world‐wide, with 32% of all global deaths in 2019. In silico studies on cardiac disease development and progression provide a tool to identify predominant mechanisms of pathological changes within the multiscale and multiphysics cardiovascular system. , Thus, in silico models are a great opportunity to identify bio‐markers and to predict disease progression and treatment outcome. , , While machine learning models can combine multiscale data from large cohorts and might identify pathological phenotypes, , , bio‐mechanical modelling provides a rich tool to explore hypotheses obtained from statistical analysis and gain understanding about underlying mechanisms.

Cardiac microstructure is a complex system of branching and inter‐connecting myocyte aggregates , assembled in myolaminae. , , Abstract representations of this complex structure, such as the fiber‐sheet model on a tissue scale, are essential to incorporate the information into models. The relation between cardiac microstructure and cardiac function renders the patient‐specific microstructure a vital component of a personalized model: (1) faster electrical conduction in the myocyte direction influences electrophysiology, (2) active stress development mainly in the myocyte direction determines cardiac contraction, and (3) anisotropic stiffness induced by the microstructure affects the passive material response. Due to this connection of cardiac electro‐mechanics and microstructure, statistical studies on cardiac function based on large data sets might profit from including representations of microstructure.

In bio‐mechanical modelling, sensitivity studies on the passive mechanical response , , , and the active contraction , showed significant influences of the underlying microstructure on stress and strains. The work by Barbarotta and Bovendeerd further suggests that models with personalized geometries but generic aggregated myocyte representations may suffer from high uncertainties in end‐systolic strains. The influence of local variations of average cardiomyocyte long‐axis orientation (ACO) compared with homogeneous structures highlights the benefit of more detailed, personalized representations of ACO. , Campos et al suggest a higher sensitivity of cardiac models to the ventricular geometry and ACO than to correct assignment of constitutive parameters. These findings demonstrate the need for patient‐specific microstructure in personalized bio‐mechanical models.

In vivo cardiac diffusion tensor imaging (cDTI) provides a non‐destructive way to measure ACO. ACO is inferred from the first eigenvector of the measured diffusion tensor. However, due to long scan times, volumetric, high‐quality measurements are challenging in clinical settings. Limited scan time requires a trade‐off between signal to noise ratio (SNR), in‐plane resolution, and the number of short‐axis imaging slices. To overcome these limitations, data interpolation strategies allow for obtaining information about ACO on a dense mesh within the left‐ventricular (LV) volume. In vivo data can be either interpolated directly or used to fit the parameters of rule‐based , , , or low‐rank models derived from observation from histology data or ex vivo cDTI data.

Rule‐based methods approximate the ACO by a function of transmural position with boundary values defined at the endo‐ and epicardial walls. , Recent rule‐based methods are based on linear functions of helix and transverse angle and use the solution of the Laplace equation with Dirichlet boundary conditions at the endo‐ and epicardial walls to determine the shape‐related transmural coordinates. , , ,

In our previous work, we proposed two low‐rank models, derived from ex vivo cDTI data, using order reduction techniques. Those models have a higher degree of freedom and are more flexible to adapt to data than classical rule‐based approaches.

Tensor interpolation strategies , provide an alternative approach to interpolate and extrapolate the tensor field to derive ACO directly from sparsely distributed data, without a model prior obtained from other data sets.

In the present work, we compare five interpolation techniques—two tensor interpolation approaches, one rule‐based approximation, and two data‐driven, low‐rank models—on five ex vivo and in vivo porcine data sets. Further, we investigate the effect of in‐plane resolution, SNR, and number of acquired short‐axis imaging slices of the cDTI input data on the interpolation result. This comparison can be used to guide future studies based on in vivo cDTI data.

MATERIALS AND METHODS

We performed two comparisons based on ex vivo data (ex vivo comparison, Section 2.3; synthetic in vivo comparison, Section 2.4) and one in vivo comparison (Section 2.5) to compare five interpolation methods from three categories: two low‐rank models, two shape‐adapted tensor interpolation approaches using two different coordinate systems, and one rule‐based method. All comparisons were performed using data from five hearts. Table 1 summarizes the three performed comparisons outlined in Sections 2.3, 2.4, and 2.5. In an ex vivo comparison, we investigate the effect of different imaging parameters of the input data. Figure 1 depicts the general workflow for all comparisons. The respective cDTI input data are pre‐processed, including segmentation, calculation of the cDTI tensors, and LV geometry generation. This is followed by mapping of the data onto a shape‐adapted coordinate system, either heat flux coordinates (HFCs) or prolate spheroidal coordinates (PSCs). These data are then interpolated to a target grid corresponding to the voxel positions, where measured ground truth is known. We interpolate and evaluate the predominant long‐axis orientation of cardiomyocytes averaged across a voxel, referred to as ACO. It coincides with the first eigenvector of the diffusion tensor. ,

TABLE 1

Summary of input and target data properties of the ex vivo, synthetic in vivo, and ex vivo comparisons

Data properties	Ex vivo comparison	Synthetic in vivo	In vivo comparison
		comparison
Number of hearts	5	5	5
Number of input data	5 (SNR variation)	1	1
settings	×2 (resolution settings)
	×7 (number of short‐axis
	input slices) = 70
Input data properties
SNR values	10, 15, 20, 25, 30,	30 (for 4 hearts)	30 (for 4 hearts)
	(was varied within the	25 (for one heart)	25 (for one heart)
	comparison for all hearts,
	respectively)
Data resolution	1.5 mm × 1.5 mm × 8 mm	1.5 mm × 1.5 mm × 8 mm	2.0 mm × 2.0 mm × 8 mm
	or		zero‐filled to
	2.5 mm × 2.5 mm × 8 mm		1.3 mm × 1.3 mm × 8 mm
Number of short‐axis	3 to 9	7 (for 3 hearts)	7 (for 3 hearts)
slices	(was varied within the	5 (for one heart)	5 (for one heart)
	comparison for all hearts,	4 (for one heart)	4 (for one heart)
	respectively)
Ground truth/target data properties
location	3D volumetric data within	one mid‐ventricular short‐	one mid‐ventricular short‐
	the myocardial mask on an	axis slice	axis slice
	isotopic grid with 0.5 mm
	resolution
SNR range	27.8 ± 7.6	30 (for 4 hearts)	30 (for 4 hearts)
		25 (for one heart)	25 (for one heart)
Data resolution	0.75 mm3 isotropic	1.5 mm × 1.5 mm × 8 mm	2.0 mm × 2.0 mm × 8 mm
	zero‐filled to		zero‐filled to
	0.5 mm3 isotropic		1.3 mm × 1.3 mm × 8 mm

FIGURE 1

Workflow. First, the cDTI input data are prepared in (I) a data pre‐processing step including (IA) data segmentation, followed by (IB) a mapping to a shape‐adapted coordinate system—HFCs or PSCs. Then, (II) interpolation of the ACO is performed, comparing five methods categorized into three interpolation strategies: shape‐based tensor interpolation, fitting a rule‐based method, and fitting a low‐rank model

Data acquisition

All data were retrospectively collected from previous studies. , Five healthy porcine hearts were imaged with a 32‐channel cardiac coil (in vivo) and an 8‐channel head coil (ex vivo), on a clinical 1.5 T MR system (Achieva, Philips Healthcare, Best, The Netherlands) with a gradient system delivering 80 mT/m per physical axis at a slew rate of 100 T/m/s. Experimental procedures were approved by the Cantonal Veterinary Office (Zurich, Switzerland) under licenses ZH072/16 and ZH 152/2013.

In vivo imaging was performed during ventilated breathing without gating, triggered to 65% of peak contraction. , Porcine anatomy differs from that of humans. The heart is not in contact with the diaphragm and hence its displacement does not follow the respiratory motion as found in humans. A second‐order motion‐compensated spin‐echo sequence , was used. Three encoding directions at  s/mm2, 3 encoding directions at  s/mm2, and 12 directions at  s/mm2 were acquired with R‐R intervals,  ms, and eight signal averages. The in‐plane resolution was 2.0 mm × 2.0 mm, zero‐filled to 1.3 mm × 1.3 mm; slice thickness was 8 mm. Two scans, with packages of five slices each, with a slice gap of 8 mm, were acquired in interleaved fashion. This resulted in a total of 10 consecutive acquired slices (no slice gap). The field of view was limited in the phase encoding direction using non‐coplanar excitation. The total acquisition duration including preparation phases was approximately 28 min (at 60 bpm).

Prior to ex vivo imaging, the hearts were arrested with a potassium chloride overdose and explanted. To ensure that the shape was preserved during fixation, both ventricular chambers were kept under hydrostatic pressure by clamping the aorta, pulmonary veins, and venae cavae. Fixation by retrograde perfusion with a 4% formalin solution and a 150 to 180 cm water column was conducted approximately 15 min after cardiac arrest. Ex vivo scans were performed with a 3D multishot diffusion‐weighted spin‐echo sequence and echo‐planar imaging readout. One unweighted (  s/mm2) image volume, 8 diffusion encoding directions at  s/mm2, and 24 directions at  s/mm2 were acquired with  s,  ms, 25 EPI interleaved profiles per shot, bandwidth 664/21 Hz (frequency/phase encoding), field of view 100 ×100× 62 mm3 to 140 ×108 × 75 mm2 depending on the size of the hearts, and two signal averages. The encoding directions for each b‐value set were computed by solving the electrostatic repulsion problem for each shell individually. The spatial resolution was 0.75 × 0.75 × 0.75 mm3, zero‐filled to 0.5 × 0.5 × 0.5 mm3. The resulting scan duration was 20 to 24 h. The ex vivo data are available at https://doi.org/10.3929/ethz-b-000478084. Cases H02, H04, H05, H06, H07, H08, PV02, and H09_2 were used to extract the low‐rank models. Cases P03, PV01, P14, P15, and P16 were used for comparison of the interpolation methods.

For both in vivo and ex vivo imaging, the system's pre‐emphasis and the product eddy‐current compensation were used to reduce eddy currents. To compensate for residual errors, non‐rigid image registration was performed to align the acquisitions with different diffusion encoding directions and b‐values.

Data pre‐processing

Data pre‐processing consisted of segmentation, generation of the LV geometry, and mapping of the cDTI data onto a shape‐adapted coordinate system.

Masks of the LV myocardium were segmented for the ex vivo and in vivo data sets. The high‐resolution, ex vivo data sets were segmented in a semi‐automatic procedure as previously reported, implemented in MATLAB 2020a (MathWorks, Natick, MA) followed by manual correction. Starting from manually selected seed points, graph‐based segmentation , was used to segment the LV lumen. To exclude the papillary muscles and trabeculae, the orientation of the first eigenvector of the diffusion tensor was used. Maps of the angle between the first eigenvector and the short‐axis image plane were determined as an approximation of the helix angle. Areas with an angle steeper than 55° and connected to the lumen were excluded from the mask of the myocardium. Edge detection was applied to obtain the endocardial boundary. The subepicardium was segmented based on a fractional anisotropy (FA) map, excluding areas with FA < 0.2 or helix angle approximation more than 70°. Edge detection and ray tracing starting at the subendocardium were combined to delineate the epicardial boundary. Manual corrections were applied to overcome segmentation errors due fixed thresholds of FA and the approximate helix angle. The in vivo data sets were manually segmented; only the short‐axis imaging slices between the apex and mitral valve were selected from the 10 acquired slices.

The LV geometry was generated from endo‐ and epicardial points extracted from the masks of the short‐axis slices for each input data set, using a marching cube algorithm from the MeVisLab toolbox (MeVis Medical Solutions, Bremen, Germany).

A physiological, shape‐adapted coordinate system, with transmural, longitudinal, and circumferential coordinates, was calculated for each LV geometry. We computed two coordinate systems, the HFC and PSC, used within different interpolation approaches: the HFC was obtained by solving a linear heat transfer problem with Dirichlet boundary conditions for each direction. , , The boundaries for the transmural coordinate were the endo‐ and epicardial walls. For the longitudinal coordinate, boundary conditions were prescribed at the LV base and a cylinder at the apex. For the circumferential coordinate, two heat transfer problems were solved, defining Dirichlet boundaries in a plane through the anterior LV‐RV intersection and the apex, and prescribing an insulation at the apex. The local flux vectors were orthogonalized and determined the local coordinate axis. The global coordinates resulted from tracking streamlines of the flux field and normalization by the total streamline length. This step compensated for the non‐linearity of the solution of the heat transfer problem, similarly to the normalization used by Paun et al and Bayer et al. The PSC is given by three angles describing the position with respect to a truncated prolate ellipsoid. The first coordinate (Θ transmural) represents the transmural depth, the second (Θ longitudinal) is the longitudinal position, and the third (Θ circumferential) represents the angular circumferential coordinate. To represent the LV in this coordinate system, an ellipsoid was manually aligned and adapted in size to best match the LV shape. A mapping of the LV geometry to the ellipsoid was obtained by registration using a symmetric log‐domain diffeomorphic demons algorithm. , ,

Ex vivo comparison

The original high‐resolution, 3D isotropic, ex vivo data sets provided the ground truth for this comparison. The input data with reduced resolution, SNR, and number of short‐axis slices were generated from these high‐quality data.

Two resolution settings were chosen to cover the range of commonly used resolution of in vivo scans: , , , , a fine setting of 1.5 × 1.5 × 8 mm3 and a coarse setting of 2.5 × 2.5 × 8 mm3. These two input data sets were obtained by down‐sampling the high‐resolution data. To this end, the data were cropped in k‐space and averaged in the through‐slice direction. During averaging, data from 16 high‐resolution slices with 0.5 mm slice thickness (seven slices basal and eight slices apical of the target slice position) were averaged to obtain one slice with 8 mm slice thickness. This was done for each slice position of the original high‐resolution data. A minimal slice gap of 8 mm was ensured while sub‐sampling slices from these data, to prevent information from the same location being included twice. The masks were resized accordingly using nearest‐neighbor interpolation.

To vary the SNR, Gaussian distributed noise was added to the down‐sampled data. The SNR within the masked myocardium was determined by dividing the mean signal of the image by the noise level. , The noise level was computed using the standard deviation within the myocardium of the difference between the two signal averages and divided by square root of two, to compensate for signal averaging. The noise was applied to the averaged data. Input data sets with SNR values of 10, 15, 20, 25, and 30 were generated.

The number of slices was reduced by equidistantly sampling short‐axis slices, covering the full range between apex and base. Three to nine short‐axis slices were used as input settings. The number of nine slices corresponds to full volumetric coverage with 8 mm thick slices of the shortest heart.

The diffusion tensors and their first eigenvectors were computed and provided the input data.

For each generated input data set, the data were interpolated to the voxel centers of the high‐resolution imaging data, applying the five methods as described in Section 2.6. The resulting ACO was compared with the first eigenvector of the diffusion tensor from the high‐resolution data set for each voxel.

Synthetic in vivo comparison

To compare the interpolation of ex vivo and the in vivo data, the in vivo comparison was replicated using down‐sampled ex vivo data, as detailed in Section 2.3. For each heart, the data set matching the SNR value and the number of short‐axis slices of the in vivo measurement were selected from the down‐sampled ex vivo data, with 1.5 × 1.5 × 8 mm3 resolution. In accordance with the procedure of the in vivo comparison, one mid‐ventricular short‐axis slice was excluded from the input data and used as ground truth reference data. Interpolation of the ACO to the voxel center of this ground truth slice was performed to compare the five interpolation methods.

In vivo comparison

For each heart, a mid‐ventricular short‐axis slice was excluded from the input data and used as ground truth reference for the interpolation comparison. The ACO was interpolated onto the voxel centers of this excluded short‐axis slice.

Interpolation methods

Tensor interpolation

We compare two tensor interpolation approaches based on the work of Toussaint et al, using the same tensor interpolation method but different coordinate systems, as detailed in Section 2.2. The interpolation using prolate spheroidal coordinates (IPSC), as presented in the original work, is compared with the combination of tensor interpolation and heat flux coordinates (IHFC).

Prior to interpolation, the DTI tensor input data are transformed into the underlying coordinates. The approximate DTI tensor at each target position is then calculated as a weighted mean of the input tensors in the log‐Euclidean sense. An anisotropic Gaussian kernel function with pre‐defined weighting matrix determines the weights w based on the distance between the target position x and the locations of the input data x :

The approximation of the DTI tensor D approx,  at a target position x is given by with N being the number of input data points. The resulting estimated ACO coincides with the first eigenvector of the approximated diffusion tensor. It is transformed back to Cartesian coordinates for comparison with the ground truth orientation. The anisotropic matrix H depends on the resolution, the SNR, and the number of short‐axis slices of the input data. The least squares optimization approach presented in the original work was applied with the PSC. It was optimized for one ex vivo heart, for both resolution setting at 14 sampling points within the parameter space given by SNR ∈[10, 30] and the number of short‐axis input slices ∈[3, 9]. To reduce the computational cost, a 2D second‐order polynomial was fit through the samples, for each spatial direction, providing the diagonal components of H for all imaging parameter settings. For the interpolation with the HFC, the components of H were scaled to compensate for the difference in the coordinate ranges by the following scaling factors: in the transmural direction, 1/(π/2), in the longitudinal direction, and 1/2π in the circumferential direction. The same anisotropy matrix was applied to all hearts.

Rule‐based method

The rule‐based method is based on the work of Bayer et al and uses the HFC. It was reduced to the LV only, not taking into account differences in the septum. The local ACO is given by two linear functions of transmural coordinate (t), prescribing the local helix angle and transverse angle

The parameters α endo, α epi, β endo, and β epi were adapted to the DTI input data by a least squares fit of the two linear functions to the helix and transverse angles of the measured ACO across all voxels.

Data‐driven low‐rank models

The two data‐driven models, proper generalized decomposition (PGD) and proper orthogonal decomposition (POD), introduced in our previous work, were extracted from ex vivo, high‐resolution DTI data sets of eight healthy hearts. The data were acquired with the same scan protocols as described in Section 2.1. This data set and the test data set used in the interpolation study are disjunct. Both models were extracted by applying order‐reduction techniques—PGD or POD—to the high‐resolution DTI data. The local ACO was projected onto the local coordinate axis of the HFC system. For each projection one low‐rank model was obtained. The final approximation of ACO is assembled combining the three models and the local coordinate axes at the target location.

To obtain the PGD model, a PGD basis was extracted from each heart separately. The PGD is based on the approach of separation of variables, thus each PGD basis function is a product of three 1D functions of transmural, circumferential, and longitudinal coordinates, respectively. A singular value decomposition (SVD) of each 1D function across hearts in a training data set provides the SVD basis functions that are combined with a set of weights resulting in the PGD model. The basis was truncated at six PGD modes, each represented by six SVD modes. When using the PGD model for interpolation, the weights were adapted to the sparse input data. To this end, an adapted PGD was applied, directly using the SVD basis functions instead of a Galerkin basis. For fewer than six input slices, the number of SVD modes for the 1D functions of longitudinal position was reduced to the number of input slices to prevent the problem from being under‐determined.

To extract the POD model, the data were first mapped onto a common, equidistant grid in HFC with 20 transmural, 200 circumferential, and 120 longitudinal points. A POD across all training hearts and transmural slices was performed and the resulting basis was truncated at eight modes. These modes were combined with a set of weights in the POD model. To adapt the weights to the sparse input data, the input data were mapped onto the same common grid. Only the grid points closest to each input data point were selected and a trilinear interpolation onto the selected grid points was performed. A gappy POD was applied to obtain the weights of the POD model.

Error metric

At each voxel center of the target data, the angular difference between the interpolated and the ground truth ACO was calculated, taking into account the sign invariance of the first eigenvector of the diffusion tensor. Due to its non‐Gaussian error distribution, the median across the target points was used for quantitative analysis.

For the ex vivo comparison, the high‐resolution (0.5 × 0.5 × 0.5 mm3), ex vivo data provided the ground truth for each heart, respectively. The target points were all voxel centers within the 3D LV mask. For the synthetic in vivo comparison and the in vivo comparison, one mid‐ventricular short‐axis imaging slice was excluded from each input data set and used as ground truth. The target points were located within this short‐axis imaging plane.

To estimate the uncertainty of the DTI measurements, the angle of the cone of uncertainty (CoU angle) , was calculated, as detailed in the appendix, from the in vivo data and the synthetic data set with same imaging parameters generated from the ex vivo, high‐resolution data. The bootstrapping required two independent acquisitions obtained by splitting the signal averages into two sets. The resulting CoU provides an estimate for half the number of averages. This reduction of the number of averages is compensated by a correction factor obtained from a simulated DTI experiment, outlined in the appendix.

To evaluate the interpolated ACO field, error maps of helix and transverse angle were determined for the in vivo comparison. The helix angle was defined with respect to the HFC. It is given by the angle between the circumferential axis and the projection of the first eigenvector of the diffusion tensor, or the ACO, to the circumferential‐longitudinal plane. The transverse angle is the angle between the projection of the first eigenvector of the diffusion tensor, or ACO, onto the transmural‐circumferential plane and the circumferential axis.

RESULTS

In the ex vivo experiment, we compared the performance of five interpolation methods for input data with different imaging parameters. Figure 2 shows the influence of in‐plane resolution, SNR, and number of short‐axis slices on the median interpolation error within the myocardium averaged across the cases. For all investigated input data settings, the interpolation error is the smallest for the two direct tensor interpolation methods (min./mean error across settings: IPSC, 10.3°/13.7°; IHFC, 9.5°/13.1°), followed by the low‐rank models (min./mean error across settings: PGD, 12.5°/18.0°; POD, 14.8°/17.8°). The RBM shows the highest interpolation error on average (min./mean error across settings: RBM, 23.7°/24.0°) with only marginal dependence on the input data parameters, as shown by a maximal error variation of 0.8° across the settings. For only three input slices and 2.5 mm 2 in‐plane resolution, the PGD model, being the most sensitive to the number of slices in the input data, results in the highest interpolation error. The sensitivity to all examined input data parameters is the highest for the PGD model and the lowest for the RBM, indicated by the standard deviation across the settings: IPSC, 2.3°; IHFC, 2.3°; RBM, 0.2°; PGD, 4.2°; POD, 2.5°.

FIGURE 2

Results of the ex vivo comparison: median angular interpolation error compared with the ex vivo, high‐resolution ground truth within the myocardium, averaged over the cases. Each subplot shows the error dependence on the number of ex vivo, short‐axis slices and the SNR of the input data. The rows present the results for the in‐plane resolution of 1.5 mm and 2.5 mm of the input data. The columns depict the five interpolation methods

Supplementary Figure 3 depicts the influence of each parameter of the input data separately as well as the standard deviation across hearts and the error histograms. The standard deviation of the interpolation error across the cases is the following: IPSC, 2.8°; IHFC, 2.2°; RBM, 2.9°; PGD, 3.0°; POD, 2.9°. It increases with coarser in‐plane resolution for IPSC, IHFC, PGD, and POD (Supplementary Figure 3a). The error distributions are shown in Supplementary Figure 3b. For both direct interpolation methods a bimodal distribution is present, reflecting the two image resolutions. The RBM shows the largest error with the smallest variation.

Changing the in‐plane resolution from 1.5 mm to 2.5 mm results in an error increase (Supplementary Figure 3c). The difference of the interpolation error between those resolutions averaged across all settings is significantly smaller for the RBM (0.3°) compared with the other methods (IPSC, 3.7°; weIHFC, 3.4°; PGD, 4.3°; POD, 3.3°). For the low‐rank models, the impact of the in‐plane resolution increases with decreasing number of input slices for fewer than six input slices.

The impact of the SNR, in the investigated range between 10 and 30, is smaller than the observed influence of in‐plane resolution and number of slices (Supplementary Figure 3d). However, the effect of the SNR on the interpolation error amplifies for SNR values below 20. The increase in interpolation error, caused by lowering the SNR from 30 to 10, is the highest for the POD model and insignificantly small for the RBM with the following error differences averaged across all settings: IPSC, 1.0°; IHFC, 0.9°; RBM, 0.1°; PGD, 1.1°; POD, 1.6°.

Similar trends can be observed when reducing the number of short‐axis input slices consecutively from nine to three. Reducing the number of slices below five for IPSC, IHFC, and POD or below seven for PGD increases the interpolation error (Supplementary Figure 3e) non‐linearly. The RBM interpolation error does not change significantly with lower number of slices. Reducing the number of slices from nine to three increases the interpolation error averaged across all settings by the following: IPSC, 3.8°; IHFC, 4.6°; RBM, 0.2°; PGD, 10.6°; POD, 5.2°.

Figure 3 shows the interpolation error mapped onto the myocardium for one exemplary heart using ex vivo input data with a resolution of 1.5 × 1.5 × 8 mm3, SNR of 30, seven short‐axis slices, concurrent with the in vivo settings used in this work. Confirming the quantitative results, the tensor interpolation methods, IPSC and IHFC, show blue regions of low errors below 12° within the majority of the myocardium, followed by the low‐rank PGD and POD models. The RBM approximation yields higher errors (>20°) in most regions of the myocardium, except for some mid‐ventricular areas. The highest errors are observed when data are extrapolated onto the apex for all methods. Further, endo‐ and epicardial surfaces show patches of elevated errors for all methods, most pronounced for the RBM.

FIGURE 3

Angular interpolation error mapped onto the myocardium for one exemplary test case. The ex vivo input data have a resolution of 1.5 × 1.5 × 8 mm3, SNR of 30 and consists of 7 short‐axis input slices, coinciding with the settings of the in vivo data used in in vivo comparison. The columns represent the results of the five interpolation methods

Average computational cost was estimated as the time to adapt the models to one input data set and reconstruct the 3D vector field of ACO on 950 000 target points. Times were measured on an eight‐core Intel Core i7‐10700K, 3.79 GHz desktop computer, for the ex vivo data. The tensor interpolation with HFC took 5.1 min, the tensor interpolation with PSC took 3.8 min, the RBM took 1.3 min, PGD 12.7 min, and POD 20 s.

The estimated SNR of the in vivo data was 30 for four hearts and 25 for one heart. Masking the myocardium between apex and base resulted in three hearts with eight short‐axis slices, one heart with six slices, and one heart with five slices. Thus, seven, five, or four input slices were used in the in vivo comparison. The average FA value within the target slices of the in vivo data was 0.36. The averaged FA value, measured within corresponding target slices of the synthetic comparison, was 0.34 (ex vivo).

To assess the uncertainty of the cDTI measurements the CoU was calculated for the in vivo comparison and its synthetic replicate with ex vivo data. This estimation requires the splitting of the signal averages into two sets, thus it is calculated on half of the averages. On average the CoU angle estimated from the data is 12.5° for the in vivo data and 10.8° for the ex vivo data with the same SNR and coverage settings. In a DTI simulation presented in the appendix, a factor of 0.69 (in vivo)/0.70 (ex vivo) to compensate for the bias due to the reduction of the number of averages was determined (Supplementary Figure 1). Applying this compensation results in a CoU with an angle estimated from the data of 8.6° (in vivo) and 7.6° (ex vivo). In a second DTI simulation presented in the appendix, we estimated the CoU angle as a function of FA and SNR (Supplementary Figure 2). The estimated CoU angle was 11.3° (in vivo) and 6.5° (ex vivo), for the average SNR of 29 and average FA of 0.3648 (in vivo)/ 0.3446 (ex vivo). This estimate from the DTI simulation is 2.7° higher (in vivo)/1.1° lower (ex vivo) compared with the estimate from the measurement including the correction factor.

Figure 4 shows the comparison of the interpolation error between the in vivo comparison, the replication of the in vivo comparison using ex vivo data (synthetic), and the ex vivo comparison. Confirming the observations from the ex vivo comparison, both the synthetic and the in vivo comparisons show the smallest interpolation errors for the direct tensor interpolation methods, followed by the low‐rank models, and the highest interpolation error for the RBM. The means and standard deviations across the cases are the following: in vivo, IPSC, 12.9°±1.9°; IHFC, 13.7°±1.1°; RBM, 24.8°±6.0°; PGD, 17.8°±5.2°; POD, 19.0°±3.3°; synthetic, IPSC, 10.6°±2.0°; IHFC, 11.2°±2.1°; RBM, 24.0°±2.8°; PGD, 14.9°±3.1°; POD, 14.1°±2.3°; ex vivo, IPSC, 11.1°±1.9°; IHFC, 10.6°±2.2°; RBM, 23.9°±2.9°; PGD, 14.4°±3.1°; POD, 15.4°±2.7°. The error averaged across the cases is higher in the in vivo comparison compared with the synthetic comparison by the following amounts: IPSC, 2.3°; IHFC, 2.5°; RBM, 0.8°; PGD, 2.9°; POD, 4.9°. The difference of the average error between the synthetic comparison and the ex vivo comparison is smaller than 1.5° for all methods.

FIGURE 4

Comparison of median angular difference within the myocardium averaged over the cases for the five interpolation methods. The blue, green, and grey bars represent the results of the in vivo comparison (Section 2.5), the synthetic in vivo comparison (Section 2.4), being a replica of the in vivo comparison with ex vivo data, and the ex vivo comparison (Section 2.4). The setting with the same SNR and number of slices as in the in vivo data set and 1.5 × 1.5 × 8 mm3 resolution is considered for the ex vivo comparison. The symbols represent the median error of the individual test cases

The in vivo comparison resulted in an interpolated vector field of ACO in one mid‐ventricular, short‐axis slice. Figure 5 shows the interpolated ACO for one exemplary heart. The measured ground truth is shown in grey; the interpolated vector field is color‐coded by the interpolation error. The qualitative evaluation confirms the best agreement between the interpolation and the measured data for the IPSC and IHFC, with a good match of spatial variations. The low‐rank model estimations show higher errors, but restricted to local areas and the boundaries. The POD underestimates the helix angle at the endo‐ and epicardial walls. The RBM results in larger patches (red) of elevated error underestimating the helix angle, especially at the subepicardium. This qualitative evaluation suggests a higher interpolation error at the surfaces, which is more pronounced at the epicardial surface. The helix angle at the surfaces is underestimated with all methods. To further visualize the resulting error distribution of the interpolated vector fields from the in vivo comparison, Supplementary Figure 5 shows error maps of the ACO, helix, and transverse angle for the same exemplary heart as shown in Figure 5. The absolute error in helix angle is higher than the error in transverse angle for most of the data points within the mid‐ventricular slice. Confirming the observations in Figure 5, the RBM and POD interpolations suffer from high errors in ACO and a helix angle of >45° in the inferior‐septal wall. This is an area with high helix angle values in the data, as shown in the exemplary helix map in Supplementary Figure 4.

FIGURE 5

Original (in grey) and interpolated (color coded by the angular interpolation error) in vivo diffusion tensors' first eigenvector field on a mid‐ventricular, short‐axis target slice, for one exemplary test case. The in vivo input data from seven input slices and SNR = 30 were interpolated onto the target slice. The columns represent the original DTI data and the five interpolation methods

The error as a function of spatial position within the ventricle is depicted in figure 6: in blue for the in vivo comparison and in green for the synthetic comparison using ex vivo data. The error shows no systematic dependence on the circumferential coordinate.

FIGURE 6

Spatial error variation within the myocardium for the in vivo data in blue and the synthetic in vivo comparison with ex vivo data in green. The rows show the circumferential and transmural variations. The columns present the results using the interpolation with PSC and HFC, the rule‐based approximation, and the low‐rank PGD and POD model estimations. Along both directions, the data were assigned to 10 bins. Each box plot represents the case‐averaged median and 25th and 75th (box) and 5th and 95th (whisker) percentiles obtained within the bin

The error as a function of the transmural coordinate shows elevated values at the endo‐ and epicardial walls (t < 0.1 and t > 0.9). A pronounced error maximum in the epicardial region (t > 0.8) is observed for the in vivo comparison for all interpolation methods. For the synthetic data, the error increase towards the epicardial wall is less distinct for the IPSC and IHFC compared with the RBM, PGD, and POD. In the epicardial region, the error is >50% higher for the in vivo comparison compared with the synthetic comparison (difference between the in vivo and the synthetic comparisons relative to the synthetic comparison for t > 0.8: IPSC, 186%; IHFC, 173%; RBM, 52%, PGD, 131%; POD, 158%). This error difference is not present at the endocardial wall.

DISCUSSION

We have compared five methods—two direct tensor interpolation strategies, two low‐rank models, and a rule‐based approach—to interpolate cDTI, short‐axis imaging data using two comparisons based on ex vivo data and one in vivo comparison. For the interpolation and extrapolation of cDTI data, the direct tensor interpolation is superior to the low‐rank and rule‐based models considered in this study. The direct tensor interpolation exploits the smoothness of the ACO field and its high spatial coherence in shape‐adapted coordinates. The comparison between the PSC and HFC shows differences in error averaged across the cases of <1.1° for all comparisons and settings. The tensor interpolation approach approximates the ACO field based on the input tensors and the coordinate system adapted to the LV shape, both specific for the individual heart. In contrast, the low‐rank models and the rule‐based method use a prior of the ACO extracted from other hearts and thus are less flexible when adapting to the input data. The rule‐based method provides only few degrees of freedom; in consequence, it cannot capture the detailed ACO field. This is confirmed by an average interpolation error of >23° and a low sensitivity of the interpolation error to the input data settings. The PGD and POD models have a higher number of degrees of freedom and can better match the target data. However, the adaption of the weights, and hence personalization to the given input, requires a minimum number of data. For the PGD model, due to the approach of separation of variables, each 1D function of the underlying basis is fit separately within each iteration of the greedy PGD algorithm. To avoid an underdetermined system of equations for fewer than six input slices, the number of degrees of freedom was reduced to the number of input slices. This truncation of the SVD basis in the longitudinal direction led to a pronounced increase in error when reducing the number of input slices below six. The POD model uses 2D modes of constant transmural depth. Consequently, the data points along the circumferential and longitudinal directions are used simultaneously to fit the weights, providing more samples per mode. Thus, the basis is not truncated for a low number of slices in the input data. However, the sensitivity of the interpolation to reducing the number of slices below five increases for the POD model and the direct tensor interpolation approaches. The larger gaps between available input points in the longitudinal direction reduce the information on the longitudinal variation given by the input data, leading to an increased error of the estimated ACO.

Areas of elevated errors for all interpolation methods are the endo‐ and epicardial surfaces and the apex. The uncertainty in the data is higher at the boundaries compared with the mid‐wall due to partial voluming and segmentation errors. A further reason for the apical error might be a potential difference in microstructure orientation relative to the LV shape compared with the mid‐myocardium. The model‐based approaches might not be able to capture the apical ACO if a higher variation between individuals exists.

In the ex vivo comparison, we investigated the influence of SNR, in‐plane resolution, and number of slices of input data on the interpolation error evaluated at the voxel positions of high‐resolution ground truth data. All these parameters can be traded off against acquisition time and are interlinked. Thus, the question arises of which parameter combinations result in similar scan times and thus can be compared to evaluate the effect of each parameter on the interpolation error. Both input data sets, with 1.5 × 1.5 mm2 and 2.5 × 2.5 mm2 in‐plane resolution, were obtained from retrospectively down‐sampled high‐resolution data. Hence, we neglect effects on SNR induced by changes in readout trajectory, voxel size, and sampling dwell time with an increase in resolution. We can assume a proportional relation between the SNR of the data and the product between voxel volume and the square root of the number of averages. Applying this estimate to the difference between the two investigated in‐plane resolutions, the volume of a voxel changes by a factor of . Consequently, this gain in resolution causes either (a) at constant scan time, a loss in SNR by a factor of , approximately resulting in a reduction from SNR 30 to SNR 10 (30 × 0.36 = 10.8); or (b) when keeping the SNR constant by ‐fold averaging, eight times longer scan time. Alternatively, this increase in scan time could be traded to acquire eight more slices at the reduced SNR. Based on the first consideration (a), we can compare the settings (1) with in‐plane resolution of 2.5 × 2.5 mm2, SNR = 30, and the settings (2) with in‐plane resolution of 1.5 × 1.5 mm2, SNR = 10. The second results in smaller interpolation errors for all interpolation strategies on average. This might result from the smoothing of the input data during interpolation and a sufficiently high information content of the interpolation input with an SNR of 10. Thus, for interpolation, the change of in‐plane resolution from 2.5 × 2.5 mm2 to 1.5 × 1.5 mm2 leads to a higher benefit than the increase in SNR from 10 to 30. It is noted that increasing spatial resolution prolongs readout duration and reduces readout bandwidth. As a consequence, motion occurring during readout may lead to spatial inconsistency of information contained in low and high spatial frequencies. To date, shortening readout duration by means of parallel imaging or multishot acquisition schemes has allowed for performance of cDTI at spatial resolutions of 1.6 × 1.6 mm2 and 1.8 × 1.8 mm2. , Additionally, lowering the readout bandwidth leads to increased susceptibility to off‐resonance artifacts that may require dedicated treatment. , Based on the second consideration (b), the additional time needed to increase the SNR from 10 to 30 could instead be used to sample additional slices. The maximal difference in the number of slices in the comparison was 6. Comparing the evaluated input settings closest to this change—(3) three slices, SNR = 30 and (4) nine slices with SNR = 10—we observe an advantage in taking more short‐axis slices into account compared with increasing the SNR. However, the improvement per added input slice reduces when adding more than six short‐axis input slices and depends on the choice of the interpolation strategy. The influence of SNR is higher for SNR lower than 20. A combination of an increase in number of input slices and SNR would be preferable.

The results presented here could be used for experimental design when a 3D volumetric ACO field should be estimated by interpolation on a fine grid or mesh. However, in practice when increasing the in‐plane resolution, the sampling changes. When keeping the dwell time constant this leads to a prolonged acquisition time by the factor of resolution increase squared. The simplified relation of SNR and resolution holds when assuming that the prolonged echo time resulting from the longer acquisition time is shorter than T 2. Given the experimental setup in this study, increasing the spatial resolution from 2.0 × 2.0 mm2 to 1.5 × 1.5 mm2 would increase echo time by ≈20 ms. Assuming a T 2 of 60 ms (Reference ) leads to an additional signal attenuation of ≈20%, which needs to be considered. Moulin et al compensated for the prolonged echo time by partial Fourier sampling and an under‐sampling factor of 2 in conjunction with a more powerful gradient system to achieve an echo time of 61 ms (versus 81 ms) for a spatial resolution of 1.6 × 1.6 mm2. Our simulation does not consider blurring due to broadening of the point spread function by decay during echo planar imaging. The effective resolution is always reduced by k‐space truncation during sampling but also due to signal dephasing during acquisition of k‐space. Further, the approximation does not consider effects of motion sensitivity and off‐resonance resulting from the longer readout.

Inspection of the results across the comparisons confirms the observed advantage of the direct tensor interpolation, followed by the low‐rank models and the rule‐based method. The ex vivo comparison mimics the application of interpolation techniques to estimate the full ACO field on a mesh of the myocardium with finer resolution than the input data. The evaluated target data points are a combination of three categories: data points sampled densely within a lower‐resolution short‐axis input slice, data points between slices, and data points outside the covered area. The same target data points were not available for the in vivo data, because the in vivo acquisition of a volumetric, high‐resolution ground truth data set is very challenging and part of ongoing research. The synthetic in vivo comparison based on ex vivo data was performed to exclude differences in the interpolation target when comparing the interpolation of in vivo data and ex vivo data. Since the ground truth data contained imaging artifacts and processing errors from registration and segmentation, an analysis with a single target slice is more likely to be affected by local errors than an analysis with target points distributed in a 3D volume. A direct comparison of the interpolated in vivo data to the high‐resolution ex vivo ground truth was not made, because it does not directly reflect the interpolation error, but might be dominated by (1) differences in the ACO due to different time points within the cardiac cycle (in vivo, systolic configuration; ex vivo, diastolic configuration); , , and (2) a potential mapping error due to the change of the geometry between the two states and thus small differences in the shape‐adapted coordinates. The in vivo comparison resulted in higher errors than its synthetic replicate from ex vivo data (averaged difference: IPSC, 2.3°; IHFC, 2.5°; RBM, 0.8°; PGD, 2.9°; POD, 4.9°). This might originate from a lower accuracy in the in vivo data despite the same SNR and resolution, due to differences during the DTI acquisition (18 versus 32 diffusion encoding directions in vivo versus ex vivo) or remaining registration errors between the single diffusion weighted images and potential slice misalignment due to residual breathing motion and zero‐order eddy currents. , The difference in the CoU angle is small for the values estimated from the data (estimates from the data, in vivo, 8.6°/ ex vivo, 7.6°; estimates from the DTI simulation, in vivo, 11.3°; ex vivo, 6.5°). However, both estimates confirm the tendency of a higher uncertainty in the data in vivo.

As shown in the spatial error analysis over the transmural coordinate, the elevated error in the in vivo comparison is pronounced at the epicardial wall for all interpolation strategies, suggesting a systematic error. In vivo, registration errors might induce uncertainties at the boundaries, which are not present ex vivo. Compared with the subendocardium, measurement errors at the subepicardium might be elevated due to the proximity to the lungs and coronary veins inducing off‐resonance effects. Segmentation errors of the subendocardium towards the papillary muscles might induce smaller errors compared with including pericardial fluid at the subepicardium. A further error source at the subepicardium, present in vivo and ex vivo, might be subepicardial fat leading to fat shift artifacts.

The evaluated error is referred to as interpolation error, even though data extrapolation to target points outside the spatial range of the input data was performed at the basal plane and apex, and all methods did not enforce that estimated ACO remains unaltered at locations where input data are given. This helps to compensate noise in the data but might lead to over‐smoothing in more heterogeneous areas. The interpolation error was evaluated at the voxel center of the measured ground truth data. The term ground truth is used even though the data reflect the ACO on a meso‐scale but not the complex microstructure on sub‐voxel scales and include measurement noise. The data quality of input data and ground truth both influence the interpolation error. Measurement errors in the input data, such as noise, perturb the resulting interpolated ACO field. Measurement errors in the ground truth reference data result in a lower limit of the interpolation errors even if the true, error free ACO field were to be evaluated. To demonstrate the quality of the ground truth data used in the comparisons, helix and transverse angle maps for a mid‐ventricular short‐axis slice are shown in Supplementary Figure 4. To quantify the measurement induced error, the uncertainty of the data was estimated by the angle of the CoU. This residual error of 8.6° (in vivo) and 7.6° (replicate from ex vivo data) in the data provides a reference of the lower bound of the error that cannot be overcome by the choice of the interpolation method. The interpolation error for the best interpolation strategy is 12.9° (in vivo) and 10.6° (replicate from ex vivo data). Thus, the interpolation error is 4.3° (in vivo)/3.0° (replicate from ex vivo data) higher than the 95% certainty range of the cDTI data.

One application of interpolation of cDTI data is patient‐specific modelling. The microstructure influences cardiac electrophysiology and mechanics. , , , , , , Interpolation is necessary to populate a finite element mesh with a patient‐specific ACO field from cDTI data. The significant influence of a heterogeneous ACO field compared with homogeneous, rule‐based representations , supports the need for realistic representations in patient‐specific models. The direct tensor interpolation and low‐rank models can capture more localized variations, not represented with the global rule‐based method. Although errors in the range of 10° influence cardiac mechanics, , , a more realistic representation of ACO leads to an error reduction in the models and thus has a substantial effect to improve image‐based, patient‐specific cardiac modelling. The interpolation error with the rule‐based method is 11.8° (on average, in vivo) higher than for the direct tensor interpolation. Despite a potential influence of measurement induced uncertainties of approximately 8.7° and the interpolation error of 12.9° (IPSC, in vivo), the direct tensor interpolation of cDTI data provides a more realistic representation of ACO. Hence, incorporating cDTI data in personalized models by direct tensor interpolation can improve patient‐specific modelling if individual data are available.

Despite the smaller interpolation errors of the two direct interpolation methods, the parametric models might be preferred for synthetic cohort studies or as parametric input in machine learning based models. Both the low‐rank models and the rule‐based approach have the advantage of reducing the dimensionality of the data and thus providing a set of parameters representing ACO. These parameters are associated with underlying characteristics of the ACO. The predominant ACO of specific pathology might be represented by specific parameter ranges, providing synthetic input data for cohort studies. These parameters could further provide a way to incorporate microstructure into machine learning models, identifying correlate characteristics across scales from large data sets. If no cDTI data are available, parametric models might enable the inference of the microstructural parameters from data of cardiac deformation patterns in an iterative FEM based optimization. However, this increases the number of parameters and thus impairs parameter identifiability.

A further application of ACO interpolation could be the estimation of strains along the ACO on the whole heart by combining the interpolated ACO field with a fiber strain estimation method based on cDTI and DENSE imaging as proposed in Reference . Uncertainties in the ACO up to approximately 13° result in small perturbations within an error range of ±0.03. The average interpolation error with the direct tensor interpolation methods from in vivo data results in errors within the same range. This suggests the feasibility of a whole heart fiber strain analysis in combination with 3D DENSE displacement fields.

A novel bio‐marker of cardiac myocardial infarction is the local propagation angle. It compares adjacent ACO vectors. 3D vector fields of the tensors' first eigenvectors are required to perform fiber tracking. Interpolation of cDTI data might support this analysis by providing ACO on a regular grid. However, dense sampling of the area of myocardial infarction is needed.

The tensor interpolation methods directly provide the interpolated tensor field. In addition to the ACO, the sheetlet direction can be obtained from the tensors' second eigenvector. This might be an advantage compared with the low‐rank and rule‐based models for studies on pathology accompanied by altered sheet angles, such as dilated cardiomyopathy (DCM), hypertrophic cardiomyopathy (HCM), cardiac amyloidosis, and aortic stenosis. , , , ,

CONCLUSION

In this study, the errors of different interpolation strategies for estimating ACO fields hve been investigated to guide the choice of an interpolation strategy in future studies requiring full LV ACO fields. The direct tensor interpolation strategies provide a better estimate compared with model‐based approaches. The ex vivo comparison provides a guidance for experimental design, suggesting that the increase in in‐plane resolution despite the associated SNR reduction improves the ACO estimate obtained by interpolation.

NBM_4667_SupplementaryMaterial_cleanVersion.docx

Click here for additional data file.

NBM_4667_supp_figure_FA_SNR_1.png

Click here for additional data file.

NBM_4667_supp_figure_FA_SNR_2.png

Click here for additional data file.

NBM_4667_supportingFigure3.png

Click here for additional data file.

NBM_4667_SuppFig4.png

Click here for additional data file.

NBM_4667_SuppFigure5.png

Click here for additional data file.