TBS: Tensor-Based Supervoxels for Unfolding the Heart

Abstract

Investigation of the myofiber structure of the heart is desired for studies of anatomy and diseases. However, it is difficult to understand the left ventricle structure intuitively because it consists of three layers with different myofiber orientations. In this work, we propose an unfolding method for micro-focus X-ray CT ($\mu$CT) volumes of the heart. First, we explore a novel supervoxel over-segmentation technique, Tensor-Based Supervoxels (TBS), which allows us to divide the left ventricle into three layers. We utilize TBS and B-spline curves for extraction of the layers. Finally we project $\mu$CT intensities in each layer to an unfolded view. Experiments are performed using three $\mu$CT images of the left ventricle acquired from canine heart specimens. In all cases, the myofiber structure could be observed clearly in the unfolded views. This is promising for helping cardiac studies.

1 Introduction

The heart consists of many muscle fibers (myofibers) useful for producing much power to pump up the blood. When the myofiber structure of one chamber, the left ventricle (LV), becomes structurally insufficient, diseases such as cardiac insufficiency may occur. More detailed understanding is desired for many purposes, such as investigating mechanisms of diseases, anatomical study, and simulation of beating. Nevertheless, investigation of myofiber structure in three-dimensional space is difficult due to the lack of analyzing and visualization techniques. Because the LV (Fig. 1(a)) has a tube-like shape and consists of three layers having different myofiber orientations [1], it is difficult to understand intuitively how myofibers are running in each layer.

Diffusion Tensor Magnetic Resonance Imaging (DT-MRI) has been commonly utilized for analyzing the myofiber structure of LV [2]. Nowadays, micro-focus X-ray CT ($\mu$CT) is also utilized [3] due to its high resolution, contrast, and low imaging cost. However, direct visualization of the $\mu$CT volumes (Fig. 1(b) and (c)) has limited usefulness for observing each myocardial layer. For intuitive understanding of each layer of the LV, we propose an unfolding method using $\mu$CT volumes. This method visualizes a $\mu$CT volume containing the LV as a stack of images showing each layer flattened.

To unfold each layer of the LV, we perform intensity mapping from each layer to unfolded images. There are several segmentation methods of layers for cardiac MRI volumes [4], which have much simpler content than $\mu$CT volumes, but no such method has been developed specifically for $\mu$CT volumes.

The three layers on $\mu$CT volumes have similar intensities as shown in Fig. 1(b). Moreover, the surrounding tissues, such as the papillary muscle, appear on $\mu$CT volumes having similar intensities. To divide the layers, we introduce supervoxel over-segmentation. Supervoxels are typically sub-volumes that consist of voxels having similar intensities, and their shapes are constrained to be regular and compact. However, we focus on the myofiber orientations shown in a stripe-like pattern (Fig. 1(b)). For this purpose, structure tensors are useful for analyzing myofiber orientations at each voxel in the $\mu$CT volumes [3]. The tensor matrix is computed from the intensities of neighboring voxels. It represents the orientations and strength of intensity fluctuations around a voxel.

Fig. 1.

A specimen of the heart shown in (a), and its $\mu$CT volume focusing on LV. Three layers having different myofiber orientations are piled up as annotated on (b) axial slice. However, any slice types shown in (c) 3D view are not useful to observe each layer.

Therefore, we propose Tensor-Based Supervoxels (TBS) by integrating a structure tensor with a supervoxel over-segmentation algorithm, instead of using the intensity values of each voxel. As illustrated in Fig. 2, each supervoxel will have homogeneous myofiber orientations, which prevents them from being simply divided into myofibers (brighter) and extracellular matrices (darker). By utilizing TBS, we can extract the middle layer by selecting supervoxels that consist of horizontal myofibers.

The technical contribution of this paper is its proposal of TBS for obtaining supervoxels having homogeneous myofiber orientations. In this paper, we show the application of TBS to the task of unfolding the heart. Although it is commonly applied to other organs such as the colon [6], for the heart there are only few works that virtually unfold the heart into 2D space [7]. While there are several technical works [8] for analyzing heart anatomy in 3D, to the best of our knowledge, this is the first work that visualize each layer by unfolding.

Fig. 2.

Idea of TBS. (a) Conventional Simple Linear Iterative Clustering (SLIC) [5]. Each supervoxel contains voxels having similar intensities. Two voxels having low and high intensities belong to different supervoxels. (b) TBS. Each supervoxel will have similar tensors. Even if voxels have different intensities, they will belong to the same supervoxel if tensors on the voxels are similar.

Fig. 3.

Overview of proposed method. (1) Supervoxel over-segmentation. We take into account the inclination angle, which can be computed with a tensor. (2) Middle-layer extraction. We select supervoxels that have a small inclination angle and then draw a smooth surface $P_\mathcal{M}$ along $\mathcal{M}$. (3) Intensity projection. We obtain intensity values on $P_\mathcal{M}$.

2 Method

2.1 Overview

The proposed unfolding method of one $\mu$CT volume I is illustrated in Fig. 3. It consists of the following steps: first, we perform (1) supervoxel over-segmentation based on tensor computation (TBS), and for each slice we perform (2) middle-layer extraction and (3) intensity projection. The output is a stack of 2D images ordered from inside to outside of the helical heart.

2.2 Supervoxel Over-Segmentation

We perform supervoxel over-segmentation of I. Supervoxels are usually segmented as if each supervoxel had a set of voxels having similar intensity values and a compact shape. However, we are not interested in the intensity value of each voxel but in the myofiber orientations. We first explain how to estimate the myofiber orientation, and then we show how to integrate it with supervoxel over-segmentation.

Estimating myofiber orientation: Myofiber orientations at a voxel $\mathbf {x}$ can be estimated as the eigenvector corresponding to the smallest eigenvalue of the structure tensor [3]. The structure tensor for a voxel $\mathbf {x}$ in a volume is typically given by

$$\begin{aligned} T(\mathbf {x})= & {} \sum _{\mathbf {x}' \in N} G( ||\mathbf {x}-\mathbf {x}'||, \sigma ) \, \mathbf {g}(\mathbf {x}') \, \mathbf {g}^\text {T}(\mathbf {x}') , \end{aligned}$$

(1)

where N represents the set of voxels around $\mathbf {x}$, $\mathbf {x}'$ represents a voxel in N, $G(||\mathbf {x}-\mathbf {x}'||, \sigma )$ represents the Gaussian function whose distance from the center is $||\mathbf {x}-\mathbf {x}'||$, and standard deviation is $\sigma$. $\mathbf {g}(\mathbf {x})$ represents the local intensity gradient around $\mathbf {x}$. Eigenvalues $\lambda _1({\mathbf x}), \lambda _2({\mathbf x}), \lambda _3({\mathbf x}) (\lambda _1({\mathbf x}) \ge \lambda _2({\mathbf x}) \ge \lambda _3({\mathbf x}) \ge 0)$ of $T({\mathbf x})$ represent the intensity fluctuation into the corresponding eigenvectors $\mathbf {e}_1({\mathbf x}), \mathbf {e}_2({\mathbf x}), \mathbf {e}_3({\mathbf x})$, respectively. The myofiber orientation is estimated as $\mathbf {e}_3({\mathbf x})$ because each myofiber bundle is a bright line, and there is little intensity fluctuation along it.

Tensor-Based Supervoxels (TBS): Simple linear iterative clustering (SLIC) [5] is a widely used supervoxel over-segmentation method similar to K-means clustering. It classifies the entire volume into supervoxels as follows: (i) As initialization, center voxels are defined as a grid pattern with an interval of s voxels. Each supervoxel S consists of one center voxel and its $s \times s \times s$ neighboring voxels at this time. For each supervoxel S, (ii) the distance between the center voxel $\mathbf {x}$ and a search voxel $\mathbf {x}'$ in $2s \times 2s \times 2s$ voxels around $\mathbf {x}$ is computed. The distance between $\mathbf {x}$ and $\mathbf {x}'$ is given by

$$\begin{aligned} d(\mathbf {x}, \mathbf {x}')= & {} \biggl \{ \underbrace{ \frac{ ( I(\mathbf {x})-I(\mathbf {x}') )^2 }{m^2}}_{\text {intensity term}} \, \, + \underbrace{\frac{ ( \mathbf {x} - \mathbf {x}' )^2 }{s^2}}_{\text {spatial term} } \biggl \}^{ \frac{1}{2} }, \end{aligned}$$

(2)

where m is the weighting parameter of the intensity term. For each voxel $\mathbf {x}'$, (iii) if $d(\mathbf {x}, \mathbf {x}')$ is smaller than the distance from the center of the supervoxel whose $\mathbf {x}'$ is currently assigned, $\mathbf {x}'$ is merged into S. (iv) The center voxel $\mathbf {x}$ is updated as the current center. (ii)–(iv) are repeated for several times.

In Eq. (1), the intensity term allows each supervoxel to have homogeneous intensities, while the spatial term allows it to have a regular and compact shape. However, we focus on comparing the tensors rather than the differences of intensity values at each voxel. We replace the intensity term as a new term obtained by comparing the structure tensors in order to cluster voxels having homogeneous tensors as illustrated in Fig. 2. While structure tensor $T({\mathbf x})$ for a voxel ${\mathbf x}$ is given by Eq. (1), we additionally introduce the structure tensor for a supervoxel S as

$$\begin{aligned} \overline{T}(S)= & {} \sum _{\mathbf {x}' \in S} \, \mathbf {g}(\mathbf {x}') \, \mathbf {g}^\text {T}(\mathbf {x}'), \end{aligned}$$

(3)

using voxels in S having higher intensity than the threshold computed by the Otsu’s thresholding, which are in myofiber regions. Here, denotes a supervoxel. The eigenvector $\overline{\mathbf {e}_3}(S)$ of $\overline{T}(S)$ corresponding to the smallest eigenvalue represents the average myofiber orientation in S.

Since we focus on myofiber orientations as explained above, we define the distance function of TBS. It has a tensor term instead of the intensity term, and a term between two tensors as

$$\begin{aligned} d'(\mathbf {x}, \mathbf {x}', S)= & {} \biggl \{ \underbrace{ \frac{ f^2(S, \mathbf {x}') }{t^2}}_{\text {tensor term}} \, \, + \underbrace{\frac{ ( \mathbf {x} - \mathbf {x}' )^2 }{s^2}}_{\text {spatial term}} \biggl \}^{\frac{1}{2}}, \end{aligned}$$

(4)

where

$$\begin{aligned} f(S, \mathbf {x}')= & {} \frac{180}{\pi } \arccos \left\{ | \overline{\mathbf {e}_3}(S) \cdot \mathbf {e}_3(\mathbf {x}') | \right\} . \end{aligned}$$

(5)

Here, $f(S, \mathbf {x}')$ is a function for comparing the tensors of S and around $\mathbf {x}$ and t is the weighting parameter of the tensor term. We define $f(S, \mathbf {x}')$ as the angle between the eigenvectors corresponding to the smallest eigenvalues. This means $f(S, \mathbf {x}')$ is aligned with the myofiber orientations. Note that when TBS is integrated into other applications, $f(S, \mathbf {x}')$ can be re-defined for the specified purpose.

2.3 Middle-Layer Extraction

We extract a curved line on the middle layer for subsequent intensity projection. Since it is known that the myofibers in the middle layer (myocardium) run horizontally [1], we obtain a segmentation of the middle layer as a set of supervoxels $\mathcal{M}$. We perform radial searches from the geometrical center voxel $\mathbf {x}_z$ on the z-th axial slice computed by binarizing the input image I. Each search is continued on a ray directing to $\psi \in \{ 0 \le \psi < 360 \}$ having a clockwise rotation from a direction vector $(-1, 0)^\text {T}$. On each search, we select a supervoxel $\hat{S}$ that satisfies

(6)

where $\mathcal{L(\psi )}$ represents the set of supervoxels touching the ray directing to $\psi$. Finally, we fit cubic B-splines [9] through the center of the supervoxels in each axial plane in $\mathcal{M}$. Each B-spline is constraint to be periodic along $\psi$. In order to obtain a smooth surface $P_\mathcal{M}$ of the myocardial middle layer in 3D, we perform second fitting of B-splines along z-direction for each $\psi$.

2.4 Intensity Projection

Expansion of the surface $P_\mathcal{M}$ on the middle layer allows us to obtain the lines on the outer layer (epicardium), while reduction gives the inner layer (endocardium). We define the depth e as the signed Euclidean distance from the center. $e<0$ represents the outer part from $P_\mathcal{M}$, while $e>0$ represents the inner part from it.

The output of the proposed method is a stack of unfolded images having various depths. By mapping intensity at (x, y) as $U_e(x,y)$ on an unfolded image $U_e$ whose depth is e, each row corresponds to the z-axis of I, while a column corresponds to the direction from the center voxel from $\mathbf {x}_z$. For generating $U_e$, we again perform searches from $\mathbf {x}_z$ on a ray directing to $\psi \in \{ 0 \le \psi < 360 \}$ for each z-th axial slice and then find a voxel whose depth is e. The intensity of I at this voxel is projected onto $U_e(z, \psi )$.

3 Experiments and Results

Dataset: We utilize three $\mu$CT volumes (Cases A, B, and C) of the LV base specimens of canine hearts which were created under the ethical approval of Kanazawa Medical University. Their specifications are shown in Table 1. We created specimens after inducing cardiac arrest by injecting KCl (20 mEq/20 ml) into the aortic root and then harvesting the heart quickly. Then we perfused the heart with a 10% formalin-neutral buffer solution. By cutting, dehydrating with ethanol, and embedding with paraffin, we created solid specimens useful for $\mu$CT imaging. Their cutting planes, which are orthogonal to the major axis of the LV, are put on the axial slice of $\mu$CT volumes. Parameters are set as: $\sigma =16$ voxels, $t=20$, and $m=240$. s is adjusted as 1 mm with respect to image resolution.

Table 1. $\mu$CT volumes for experiments: cases A, B and C.

Fig. 4.

Supervoxels on case A. Each supervoxel is painted as one color representing the inclination angle of myofiber orientation, which is estimated as explained in 2.2. Myofiber structure is represented as stripe-like pattern with strong contrast. Conventional SLIC supervoxels shown in (b) are jagged because they respects the local intensity of each voxel. Since TBS shown in (c) are focusing on local tensors, it divided the volume by generating more spherical supervoxels.

TBS: We compared conventional SLIC and TBS by applying them to Case A. As shown in Fig. 4, the supervoxels of conventional SLIC were jagged. Although parameter setting $m=240$ for conventional SLIC usually preserves compactness for most types of dataset, because of strong contrast and noise of $\mu$CT images, compact supervoxels could not be created. Because conventional SLIC depends on the intensities of each voxel, the supervoxels were divided into brighter (myofibers) and darker (extracellular matrices) parts. In contrast, TBS did result in more compact (spherical) supervoxels, and the three layers were effectively divided by the myofiber orientations. For evaluating compactness quantitatively, we computed sphericity of SLIC and TBS on Fig. 5. The sphericity is defined as $\pi ^{1/3} (6V)^{2/3}/F$, where V and F represent the volume and surface area of a supervoxel, respectively. TBS gives much higher sphericities than conventional SLIC, which is important to achieve a smooth and stable division.

Fig. 5.

Sphericities of SLIC and TBS.

Fig. 6.

Unfolding results of case A. Starting voxels and unfolding directions are roughly indicated as arrows on (a) input volume and (c) unfolding results. (b) shows the depth from middle layer by color. In (d) we show unfolding results at several depths e.

Unfolding: Our proposed method unfolded the LV successfully for all three cases. The results of Case A are shown in Fig. 6. From the input volume (Fig. 6(a)), we extracted the layers (Fig. 6(b)) and unfolded the LV as if straightening (Fig. 6(c)). Unfolded images with $e \in \{ -50, 0, +50 \}$ as inner, middle, and outer layers of the LV, respectively, are shown in Fig. 6(d). The different myofiber orientations of each layer as known from cardiac anatomy [1] can be well observed in the unfolded views.

4 Discussion and Conclusions

The usefulness of our novel supervoxel over-segmentation method, Tensor-Based Supervoxels (TBS), was demonstrated for unfolding the LV. TBS are potentially useful for applications that do not focus on local intensity values but rather on the types of structures, e.g. region segmentation having homogeneous intensity fluctuation. When TBS is integrated into other applications, $f(S, \mathbf {x}')$ can be re-defined for the specified purpose.

This work was carried out based on the anatomical knowledge that the LV consists of three layers: endocardium, myocardium, and epicardium. Using Tensor-based Supervoxels (TBS), the three layers were effectively divided as shown in Fig. 4. The scheme of the proposed method is also promising for segmentation of each layer. Moreover, in Fig. 6, it is shown that each layer has specific tendency of myofiber orientation, as studied in cardiac anatomy [1]. A limitation can also be observed in these figures. More improvement of middle layer extraction is desired since unfolding results have some distortions along the Z-axis. Distortion may be reduced by adjusting the width of the unfolded images with respect to the circumference of the layer. Despite this limitation, we showed that virtual unfolding of the helical heart is promising for detailed investigation of myofiber structures. It could improve studies of cardiac anatomy and diseases. Also, this will be also useful for educational purposes. Future work includes the intensity projection onto the 17 segment model of the American Heart Association, which is more common way to visualize the LV. Also, another supervoxel method using both intensity and tensor information may be useful for determining myofiber and extracellular matrix regions.