Denoising Moving Heart Wall Fibers Using Cartan Frames

Abstract

Current denoising methods for diffusion weighted images can obtain high quality estimates of local fiber orientation in static structures. However, recovering reliable fiber orientation from in vivo data is considerably more difficult. To address this problem we use a geometric approach, with a spatio-temporal Cartan frame field to model spatial (within time-frame) and temporal (between time-frame) rotations within a single consistent mathematical framework. The key idea is to calculate the Cartan structural connection parameters, and then fit probability distributions to these volumetric scalar fields. Voxels with low log-likelihood with respect to these distributions signal geometrical “noise” or outliers. With experiments on both simulated (canine) moving fiber data and on an in vivo human heart sequence, we demonstrate the promise of this approach for outlier detection and denoising via inpainting.

1 Introduction

Popular methods for denoising diffusion weighted (DW) images use a non-local means filter adapted to Rician noise [3, 4] or reduce the contribution of the less significant components in a principal component analysis-based reconstruction of the DW images [8]. Both classes of methods give impressive results and are now standard in pipelines for pre-processing diffusion images prior to diffusion tensor (DT) or high angular resolution diffusion imaging (HARDI) reconstruction. Typically such denoising methods have been applied to in vivo diffusion imaging of the brain or to ex vivo diffusion imaging of fixed tissue, such as the mammalian heart. In both scenarios, the diffusion sequence allows for high enough magnetic field strength and scan times, which when combined with denoising reveals high quality fiber orientation estimates. For instance, in the DT reconstruction of an ex vivo rat heart in Fig. 1 (left), one sees a smooth coherent rotation of fibers from epicardium to endocardium in both the left and right ventricles, along with the vertical orientation of fibers in the papillary muscles.

In the context of in vivo diffusion imaging of the human heart, which is an area of active research and ongoing progress [5, 11, 12], the recovered fiber orientations may be locally incoherent, despite DW denoising being applied prior to tensor fitting, with the spatial resolution being poor, as illustrated for the left ventricle in Fig. 1 (right). In the present article, we address the problem of detecting outliers in in vivo mammalian diffusion tensor reconstructions from a geometrical perspective. We do so by introducing a 4D (spatio-temporal) Cartan frame, which, when fit to measured data, allows us to characterize the distributions of Cartan form parameters. Using these distributions, one can then determine the degree to which a particular local frame fit to the data is predicted by the estimated distributions. This allows us to both detect outliers, and then remove and in-paint the missing regions. Our experiments on simulated moving fibers in a canine heart wall from the STACOM 2014 LV mechanics challenge show that our statistical approach is capable of robustly and accurately identifying areas with incoherent fiber orientations. We also demonstrate the applicability of our algorithm to human in vivo heart wall fiber data.

In vivo cardiac diffusion imaging is a powerful means for bringing methods that were once only applicable to ex vivo excised hearts to bear on the domain of the living. However, it faces significant computational challenges, especially for an organ continuously in motion, such as the heart. By considering both spatial and temporal information for cardiac fiber denoising in a single consistent framework, our approach shows promise in addressing these challenges.

Fig. 1.

DT reconstruction of fiber orientations with tractography run from a local region. Left: Ex vivo rat heart (LV and RV). Right: In vivo human heart (LV).

2 Cartan Forms for Moving Fibers

The attachment of frame fields to diffusion MRI data of ex vivo hearts was proposed in [9] for the parametrization of the static geometry of heart wall myofibers. We build on this construction by proposing a Cartan frame for spatio-temporal (in vivo) heart wall data, as described previously [10]. We assume that for any query time $t \in \mathbb {R}$ we have a stationary state of an orthonormal moving frame field $F^t$ attached to heart wall fibers, with the corresponding universal coordinate $E^t = [e_1, e_2, e_3]^T$: $f_1$ is aligned with the fiber orientation, $f_3$ is given by the component of the normal to the heart wall that is orthogonal to $f_1$ and $f_2$ is their cross-product. Our goal is to parametrize the local rotation of $F^t$ in both space and time so we extend the existing 3D coordinate system to add a 4th dimension to represent the time axis.

Fig. 2.

The application of a specific set of 4D connection parameters $c_{ijk}$ to a frame field attached to a fiber direction.

Let $E = [e_1,e_2,e_3,e_4]^T$ represent our extended universal coordinate system, $\mathbb {R}^4$, in which $e_i\cdot e_j = \delta _{i,j}$ where $\delta _{i,j}$ is the Kronecker delta, and $e_4 = [0, 0, 0, 1]^T$ is a basis vector for the time axis in our 4D representation. Thus a query $<x,y,z,t>$, where $x,y,z,t \in \mathbb {R}$, represents a query $<x,y,z>$ in $E^t$. We extend the local frames $F^t$ to 4D as follows:

Definition 1

Let $F^t=[f_1,f_2,f_3]^T$, where $f_i = (f_{i,x}^t,f_{i,y}^t,f_{i,z}^t)$, describe the cardiac frame field in $\mathbb {R}^3$ at time t. Then: $F=[f_1,f_2,f_3,f_4]^T$, where $f_i = (f_{i,x}^t,f_{i,y}^t,f_{i,z}^t,0) \; \forall \; i\in \{1,2,3\}$ and $f_4 = (0,0,0,1)$, is the 4D extension of the frame fields $F^t$.

The next step is to extend the existing connection forms from 3D to 4D. Let $F = [f_1,f_2,f_3,f_4]^T$ be an arbitrary frame field on $\mathbb {R}^4$, which, for $i\in \{1,2,3,4\}$, has the following parametrization in the extended universal coordinate system $E = [e_1,e_2,e_3,e_4]^T$: $f_i=\alpha _{i1}e_1 + \alpha _{i2}e_2 + \alpha _{i3}e_3 + \alpha _{i4}e_4.$ A covariant derivative of this frame field with respect to a vector v at point p is given by $\nabla _vf_i=\omega _{i1}(v)f_1(p)+\omega _{i2}(v)f_2(p)+\omega _{i3}(v)f_3(p)+\omega _{i4}(v)f_4(p),$ where $[\alpha _{ij}]=[f_i\cdot e_j]$ is the attitude matrix in $\mathbb {R}^4$ and $\omega _{ij}=\sum _{k}(d\alpha _{ik})\alpha _{kj}$ is a 1-form.

Definition 2

Let $i<j$ and $i,j,k\in \{1,2,3,4\}$. By feeding the frame field’s unit vectors ($f_i$s) to $\omega _{ij}(v)$ we have $c_{ijk}=\omega _{ij}(f_k)$.

Remark: For $j = {4}$, $c_{ijk} = \omega _{ij}(f_k) = 0$. For a point p and with $i<j\in \{1,2,3\}, \omega _{ij}(f_k)$ represents the amount of $f_i(p)$’s turn toward $f_j(p)$ when taking a step towards $f_k(p)$.

Given the skew-symmetry property for the $4\times 4$ connection form matrix we build $3\times 4 = 12$ different non-zero and unique $c_{ijk}$s. We can use these coefficients to estimate the motion of the frame field using first order approximation of the Taylor expansion of the frame field F in the direction of the vector v at point p. Then, by applying Cartan’s structural equation, with $v=\sum _k v_k f_k$, we get $f_i(v) \simeq f_i(p)+ \sum _j\left( \sum _k v_k c_{ijk}\right) f_j(p).$

Our spatio-temporal Cartan frame construction illustration in Fig. 2 shows an initial fiber direction with an attached frame field, and then applies a specific set of $c_{ijk}$s to it. Here $f_1$ is in the direction of the fiber, $f_3$ is in the in-page direction orthogonal to $f_1$ and $f_2 = f_3 \times f_1$. The figure (left to right) shows three samples in time of the orientations in the local neighborhood of the fiber, generated with the parameters $c_{123} = 0.5$ radians/voxel, $c_{124} = 0.03$ radians/time-step, and $c_{ijk} = 0$ for all remaining connection parameters. The positive $c_{123}$ value results in a clockwise rotation of fibers in the in-page direction (panel A) and the positive $c_{124}$ value results in an increase in the total in-page rotation of fibers in time (panels B and C).

We now show how to calculate the 1-form coefficients of a given 4D frame. From the definition of the $c_{ijk}$s in [9] and their first order Taylor expansions, for an arbitrary vector v at point p, and $i,j,k,n \in \{1,2,3,4\}$ we can re-write the equation from Definition 2 as $c_{ijk} = \omega _{ij}(f_k) = {f_j}^T J(f_i)f_k,$ where $J(f_i) = [\frac{\partial f_{ij}}{\partial x_k}]$ is a Jacobian matrix. Given a discretized frame field, we then apply this equation to calculate the connection form coefficients.

2.1 Connection Form Distribution Fitting for Outlier Detection

Our goal is to use connection form measurements for detecting locations where the estimated fiber orientations may be viewed as outliers. To this end, we fit PDFs to the set of observed values of the forms across voxels. To simplify our calculations, we treat each connection parameter separately; i.e. if $\widehat{C}_{ijk} = \{ \hat{c}_{ijk}(\varvec{p}) \}$ is the set of observed connection form measurements for $c_{ijk}$ across positions $\varvec{p}$, we estimate a PDF $f_{ijk}(c | \theta )$ via fitting to $\widehat{C}_{ijk}$. We considered both non-parametric fitting, via kernel density estimation using Gaussian kernels for human data, and maximum likelihood fitting with a normal distribution for simulated data, based on observations in [9]. This provides 12 independent distributions, with which we can compute the probability density of a given observed connection form value $\hat{c}_{ijk}$ via $f_{ijk}(\hat{c}_{ijk} | \theta )$.

To detect outliers, we first estimate the probability distributions described above. Then, for each set of observed connection parameter values $\widehat{C}_{\varvec{p}} = \{ \hat{c}_{ijk}(\varvec{p}) \}$ at a given voxel position $\varvec{p}$, we estimate a log-likelihood of the given observed Cartan measurements at that position via $\mathcal {L}_\theta (\varvec{p}) = \sum _{ \hat{c}_{ijk} \in \widehat{C}_{\varvec{p}} } \log ( f_{ijk}(\hat{c}_{ijk} | \theta ) )$ where we have assumed independence, in keeping with the previous assumption. We may then label the fiber vector $f_1(\varvec{p})$ at a given position as an outlier based on the approximate likelihood $\mathcal {L}_\theta (\varvec{p})$ of its observed connection values.

Once the outliers are labeled, we can perform denoising by removing the noisy fibers and using an inpainting algorithm to re-estimate them. Here, we use an approach that separately fills in the x, y, and z components of the fiber vectors directly, utilizing the global 4D data of each component as a scalar field [13].

3 Experiments

We construct two data sets for evaluating our method. The first uses canine heart data from the STACOM 2014 LV mechanics challenge [2]. This includes a hexahedral mesh at the beginning of the beat cycle, the associated local cardiomyocyte fiber orientations from the first principal eigenvector of ex vivo diffusion tensor MRI data, an endocardial pressure curve, and the positions of several reference points on the base of the left ventricle at three points in the beat cycle. We used this data for a finite element simulation of the heart wall using the transversely isotropic form of the Holzapfel-Ogden constitutive equations [6], implemented in FEBio [7] as a plugin. We further apply to each fiber a controlled rotation $\theta (d,t)$ of the undeformed fiber orientation at each point in the heart wall, as a linear function of time t and transmural distance d from the midwall, about the transmural axis (with maximal rotation of 20$^\circ$).

The second is from in vivo diffusion imaging of an entire human heart. Data of a single healthy volunteer was acquired on a 3T scanner (Philips, Achieva) using a 32-channel cardiac coil. DWI was performed using a SE sequence with cardiac triggering in free breathing with asymmetric bipolar gradients [11] and additional compensation for the slice and readout gradients. Data was acquired with 150 and 220 ms delays after the cardiac trigger and b-values of 0, 10, 20, 30, 50, 100, 200, and 400 s/mm$^2$ with 6, 3, 3, 3, 3, 3, 3, and 24 gradient-directions, respectively. The imaging parameters were: TE = 62 ms; TR = 14 heart beats; FOV = 280$\,\times \,$150 mm$^2$ (using outer volume suppression using rest slabs); slices = 14; voxel size = 7$\,\times \,$2.5$\,\times \,$2.5 mm$^3$; acquisition matrix = 112$\,\times \,$48 SENSE factor = 2.5; partial Fourier = 0.85; EPI bandwidth = 42 Hz/pix; averages = 1; fat suppression = SPAIR; Gmax = 62 mT/m; max slope = 100 mT/m/ms and acquisition time = 13 min. Data processing was done using DTITools for Mathematica 10 and included registration to correct for subject motion and eddy current deformations, principal component based noise suppression and tensor calculation using weighted linear least squares estimation. The primary eigenvector of the tensors within the left ventricle were extracted and used for subsequent analysis.

3.1 Outlier Detection with Canine Simulation Data

We first consider the utility of our detection framework using the finite element simulations of canine beating heart data described above. The simulation itself is considered to be the ground-truth, to which we compare our outlier detector. Given a series of moving heart fiber volumes in time, we first add synthetic noise to the fiber field $f_1$, and then estimate the connection form parameters at each voxel. The likelihood methodology described above is then used to label statistical outliers as noise. Given that the noise is artificial in nature, we can empirically measure the detection performance, and compare it to two alternative approaches.

Artificial Vector Noise Addition. To test the capabilities of our detection algorithm, we added synthetic noise to the fiber fields of our simulated canine datasets, via two different approaches. Each approach is applied separately to the heart time slices. The first method utilized the von Mises-Fisher (VMF) distribution, which describes a probability density function (PDF) on $S^2$ in $\mathbb {R}^3$, sampled via the Ulrich-Wood rejection sampling approach. There are two parameters: $\kappa$, which inversely controls the noise level per position, and $\rho$, the probability that a given voxel will have noise added to it. Given $\kappa$ and $\rho$, for each position $\varvec{p}$, with probability $\rho$, we replace $f_1(\varvec{p})$ with a unit vector v from the VMF distribution with density $f(v|\kappa ,\mu )$, where $\kappa \ge 0$ is the concentration parameter and the mean direction is $\mu = f_1(\varvec{p})$.

The second type of noise simultaneously alters blocks of data in the fiber field, via a random rotation matrix. Given a fixed number of noise blocks $n_B$, for each block B, we generate a random rotation matrix $R_B(\theta ,\varvec{a}) \in SO(3)$, which rotates a given vector by an angle $\theta$ about an axis $\varvec{a}$. This is done by sampling $\theta$ uniformly from the fixed interval $[\theta _{\text {min}},\theta _{\text {max}}]$ and choosing $\varvec{a}$ as a uniformly random vector from $S^2$. We then replace every fiber vector $f_1(\varvec{p})$ with its rotated counterpart $R_B(\theta ,\varvec{a}) f_1(\varvec{p})$, for every $\varvec{p}\in B$.

Comparative Approaches. We compare our detection algorithm to two alternative approaches: a rule-based approach and a local deviation detector. As in the likelihood approach, each method assigns a scalar value per voxel, which can be used to rank and thus label the associated fiber vector as noise.

The first measures variation from the fiber vector field generated by the rule-based approach, using the methods for the left ventricular wall described in [1]. The $\alpha _{\text {endo}}$ and $\alpha _{\text {epi}}$ were optimized via grid search, to fit the simulated data at $t=0$ (with total average undirected angular error $\sim$18$^\circ$). In this case, a position $\varvec{p}$ is assigned ${\varTheta }_{\text {rule}}(\varvec{p}) = \delta _{\text {err}}(f_1(\varvec{p}), f_{1, \text {rule}}(\varvec{p}))$, where $\delta _{\text {err}}(u,v)$ is the undirected angular error between u and v.

The second is a local method which labels a sudden deviation of the orientation of a fiber vector from that of its neighbours as noise. Per voxel position $\varvec{p}$, we consider its spatial neighbourhood $N(\varvec{p})$ and assign to $\varvec{p}$ the average local error $\bar{\varTheta }(\varvec{p}) = |N(\varvec{p})|^{-1} \sum _{q\in N(\varvec{p})} \delta _{\text {err}}(f_1(\varvec{p}), f_1(\varvec{q}))$.

Noise Detection Results. We demonstrate that our statistical connection form likelihood approach performs better than the rule-based or local detector algorithms (see Table 1). In particular, the method is more accurate on the block based noise, which may be interpreted as a form of local structure detection. The neighbourhood approach has difficulty with this case, because the rotational noise is consistent throughout the block. For the VMF noise, the local approach slightly outperforms the statistical approach when the concentration parameter $\kappa$ is small, but the statistical approach is better at high concentrations (small deviations). In general, the rule-based approach performs better when the noise is more severe; it has the advantage of not relying on the measured orientations, but this also means it is unable to utilize local information. Qualitatively, one can see the results of outlier detection with our algorithm in Fig. 3.

Table 1. Performance comparison of detection algorithms across noise types. For VMF noise, $\rho =0.35$. For block noise, $n_B=10$ with blocksize of $5\times 5\times 5$, and we fix $\theta _{\text {max}}-\theta _{\text {min}}=\pi /4$ and define $\bar{\theta }=(\theta _{\text {max}}-\theta _{\text {min}})/2$. Values shown are detection accuracies as percentages. There are $\sim$10$^6$ fiber vectors (i.e. unmasked voxels). For the block noise, we choose a detection threshold such that $\sim$5% of the voxels are chosen; for the VMF noise, the proportion $\rho$ is chosen.

Fig. 3.

Visualization of outlier detection method on Canine data. Top row: heart slice with VMF-based noise. Bottom row: heart slice with random block rotation noise. Left column: ground truth. Middle column: outliers detected by our algorithm in red. Right column: slice with outliers detected, removed and inpainted.

Fig. 4.

Visualization of outlier detection method on in vivo human data. Left: heart slice with outliers labeled as noise by our algorithm coloured in red. Right column: slice with outliers detected, removed and inpainted.

3.2 Outlier Detection with Human in vivo Data

Our methodology is motivated in part by the significant noise level observed in human in vivo DTI, which we hope to reduce by our outlier detection methodology. Since ground-truth orientations are not known our results are qualitative in the sense that we hope to correct the most significant errors in the measurements, to obtain a more biophysically plausible reconstruction. In Fig. 4, we show a slice of the fiber vector field, with the noise detected by our statistical algorithm (left), as well as the denoised version (right), which shows an increase in smoothness and consistency.

4 Conclusion

Herein, we have devised an algorithm for outlier detection via statistical analysis of observed Cartan connection form parameter values. Its performance and applicability to denoising have been demonstrated in both a synthetic setting, via artificial additive noise in heart beat cycle simulations of canine heart data, and for in vivo human DT measurements over time. The use of 4D connection forms represents a principled method to smooth and denoise moving fiber data in a spatiotemporally consistent manner. Such an approach can complement the algorithms presently in use, particularly for the lower signal-to-noise ratio of in vivo human DTI.