Keywords

1 Introduction

Blood vessel segmentation is a very basic and important problem in medical image processing. Blood vessel regions in different tissues have different features. A variety of segmentation methods were proposed to make full use of those different features. While thresholding and region-growing approaches utilize intensity-based criteria, Hessian-based methods utilize the second-order image derivatives to represent high-order local geometric information [1, 2]. Other methods such as model-fitting-based methods [3] and Markov chain Monte Carlo approaches [4] have achieved excellent results in blood vessel segmentation.

Partial nephrectomy is a common treatment for kidney cancers. For surgical planning of partial nephrectomy, renal artery segmentation plays an important role. Renal arteries have lower contrast in comparison with other blood vessels such as hepatic arteries and carotid arteries. Precise renal artery segmentation remained a big challenging. In previous work [5], we presented a two-step segmentation method, which utilized graph-cut to segment thick blood vessels and then utilized a template model to track and extract tiny blood vessels in deeper levels of the vascular tree. Although graph-cut show high robustness in segmentation tasks, it remains difficult to segment long thin tubular structures in low contrast with 1–2 voxels in diameter. This possibly stems from the fact that the optimal graph still cuts across the high-cost regions (i.e. the blood vessels) that have lower energy cost than even a cut across low-cost regions. In our previous work, a new Hessian-based dissimilarity measure was proposed to adjust the energy function to segment tiny blood vessels. The dissimilarity of two tensors is measured utilizing an Euclidean distance metric.

In this paper, we propose a tensor-based graph-cut method in a Riemannian metric space. Hessian matrix analysis is a valuable approach to represent a local geometric structure. We compute the Hessian matrix for each voxel as a “Hessian tensor”. We use the Riemannian metric as a tensor measure to calculate tensors statistics. This geometric information will then be used in a Gaussian Mixture Model (GMM) to estimate the probability distribution of foreground and background regions. The main contributions of this work can be summarized as: (1) accuracy improvement of graph-cut method for small blood vessel segmentation by means of tensor analysis and (2) applying Riemannian metric to calculate the Hessian tensor statistics.

2 Method

2.1 Generation of Tensor Field

Han et al. [6] proposed a multi-scale nonlinear structure tensor (MSNST) space to extract texture features from color images. Firstly, a multi-scale structure tensor of form \(\mathbf T _s = \alpha ^{-2s}\sum ^{H}_{h=1}{} \mathbf T _h\) must be constructed, where \(\mathbf {T}_h\) is the structure tensor of the h-th color channel and H is the total number of color channels. s denotes the scale. Here, \(\alpha \) is a positive constant that can be set as 2. Then, nonlinear diffusion is applied to all scales of structure tensor \(\mathbf T _s\), and finally a MSNST space is obtained as \(\mathbf \Gamma = \{ \hat{\mathbf{T }}_0, \hat{\mathbf{T }}_1,..., \hat{\mathbf{T }}_{s-1} \}\).

Instead of using structure tensor, we utilize the Hessian matrix to generate a tensor representation at each voxel. This is because Hessian matrix is more appropriate to represent the blood vessel structures. Different from texture features, the most appropriate scale for each voxel to represent vessel structures can be obtained by the optimal response of Hessian-based vesselness filter [1]. The Hessian matrix is given as

$$\begin{aligned} \nabla ^2 I(\mathbf x ) = \begin{bmatrix} I_{xx}(\mathbf x )&I_{xy}(\mathbf x )&I_{xz}(\mathbf x ) \\ I_{yx}(\mathbf x )&I_{yy}(\mathbf x )&I_{yz}(\mathbf x ) \\ I_{zx}(\mathbf x )&I_{zy}(\mathbf x )&I_{zz}(\mathbf x ) \\ \end{bmatrix} \end{aligned}$$
(1)

which can be regarded as a second-order tensor \(\mathbf {T}_H\). Here, \(I_{ij} = \dfrac{\partial ^2}{\partial i\partial j}I, (i,j = \{x, y, z\})\) represents the second-order partial derivatives of the local image \(I(\mathbf {x})\).

In contrast to the MSNST scheme, a multi-scale Hessian-based vesselness enhancement filter is utilized here to find the most appropriate scale for each voxel. Calculating the eigenvalues of the Hessian matrix, a vesselness measure can be obtained by determining the largest response among all scales. We use Sato’s vesselness measure \(\mathcal {V}\) [1] in this paper. Due to the fact that the scale with the highest response \(\mathcal {V}_{max}\) is the best scale for local representation of tubular structures at each voxel, the tensor is calculated at the scale with the highest vesselness response.

Although, every voxel of image \(I(\mathbf {x})\) has a scale at which it achieves its largest response of vesselness filter, most voxels actually do not belong to a vessel. Therefore, it is particularly important to reduce these “noisy tensors”. In this work, we replace a diagonal tensor \(\mathbf {T}_D\) with the tensor of the non-vessel voxels that can be expressed as

$$\begin{aligned} \begin{aligned} \mathbf {T} = {\left\{ \begin{array}{ll} \mathbf {T}_H, \mathcal {V} > 0 \\ \mathbf {T}_D, \mathcal {V} \le 0 \\ \end{array}\right. } \end{aligned} \end{aligned}$$
(2)

where \(\mathbf {T}_D = \begin{bmatrix} \lambda _1&0&0 \\ 0&\lambda _2&0 \\ 0&0&\lambda _3 \\ \end{bmatrix}\) and \(\lambda _3 \gg \lambda _2 \approx \lambda _1 > 0\). \(\mathbf {T}_D\) behaves as a plate-like structure in vesselness terms. An example of tensor field is shown in Fig. 1, where only tensors with \(\mathcal {V} > 0\) are illustrated. The tensor visualization technique was presented by Barmpoutis et al. [7].

2.2 Tensor Metric in a Riemannian Space

Numerous studies in the diffusion tensor imaging (DTI) literature have developed methods to perform statistics over tensor space [8, 9]. A key point of tensor computing is that the tensor space does not form any vector spaces, thus standard linear statistical techniques cannot be applied. The space of tensors is a type of manifold referred to as Riemannian manifold. The Riemannian manifold \((\mathcal {M},g)\) consists of linear subspaces of Euclidean space, the geodesic distance on a Riemannian manifold is a continuous collection of inner products on all tangent spaces \(T_x\mathcal {M}\) over the manifold \(\mathcal {M}\), known as Riemannian metric. The Riemannian metric g makes it possible to define geometric notions on Riemannnian manifolds such as the geodesic distance, the mean tensor value, and the gradient descent.

In local coordinates, a Riemannian metric is symmetric positive definite (SPD) matrix. Therefore, as the Hessian matrix may take other forms than positive definite (not positive definite, NPD), it cannot be mapped onto a Riemannian manifold. Here, we transform the NPD Hessian matrix to a SPD matrices, meanwhile preserving the structure feature. Suppose that the tensor \(\mathbf {T}^-\) is a NPD matrix. Further suppose that \(\mathbf {U}\) is an invertible orthogonal matrix with columns correspond to eigenvectors. Then \(\mathbf {T}^- = \mathbf {U}\mathbf {D}\mathbf {U}^{-1} = \mathbf {U}\mathbf {D}\mathbf {U}^T\) is diagonalization of the matrix \(\mathbf {T}^-\). Here, \(\mathbf {D} = \text {diag}(d_i)\), where \(d_i\) is the i-th eigenvalue.

Taking the absolute value of \(\mathbf {T}^-\), we will obtain a tensor \(\mathbf {T}^+\) in space of SPD matrices \(Sym^+\) as follows:

$$\begin{aligned} \mathbf {T}^+ = abs(\mathbf {T}^-) = \mathbf {U}\text {diag}(abs(d_i))\mathbf {U}^T, \mathbf {T}^{+} \in Sym^+. \end{aligned}$$
(3)

Since the geometric structure is related to the ratio of absolute value of the eigenvalues, the absolute value does not alter the geometric information, i.e. it does not alter the geodesic distances. Notice that the original positive definite Hessian matrices belong to dark structures such as dark vessels (\(\lambda _3\ge \lambda _2>\lambda _1>0\)) and dark blob-like structures (\(\lambda _3\approx \lambda _2\approx \lambda _1>0\)). These matrices have already been replaced by diagonal tensors, hence the absolute value operator does not affect the segmentation result.

According to the definition in [8, 9], the distance between two tensors \(\mathbf {T}_1\) and \(\mathbf {T}_2\) is given by

$$\begin{aligned} d(\mathbf {T}_1, \mathbf {T}_2) = \left( \sum ^{n}_{i=1}\log ^2\lambda _i(\mathbf {T}_1, \mathbf {T}_2) \right) ^{\tfrac{1}{2}}, \mathbf {T}_1,\mathbf {T}_2 \in Sym^+ \end{aligned}$$
(4)

where \(\lambda (\mathbf {T}_1, \mathbf {T}_2)\) returns the generalized eigenvalue of the tensors \(\mathbf {T}_1\) and \(\mathbf {T}_2\). and n denotes dimension which is 3 here.

The Fréchet mean is a generalization of the mean values in an arbitrary metric space that is established by minimizing the sum of squared distances: \(\bar{\mathbf {T}} = \arg \min {\dfrac{1}{N}}\sum ^N_{i=1}dist(\mathbf {T}_i, \bar{\mathbf {T}})\). The mean value in a Riemannian metric space is uniquely defined [10]. A Newton gradient descent algorithm is used to solve this minimization problem. Initially, set \(\bar{\mathbf {T}}_0 = \mathbf {T}_0\). The mean tensor at step \(t+1\) is given by

$$\begin{aligned} \bar{\mathbf {T}}_{t+1} = \exp _{\bar{\mathbf {T}}}\left( \frac{1}{N}\sum ^N_{i=1}\log _{\bar{\mathbf {T}}_{t}}(\mathbf {T}_i) \right) = \bar{\mathbf {T}}^{\tfrac{1}{2}}_t\exp \left( \frac{1}{N}\sum ^N_{i=1}\log (\bar{\mathbf {T}}^{-\tfrac{1}{2}}_t\mathbf {T}_i\bar{\mathbf {T}}^{-\tfrac{1}{2}}_t) \right) \bar{\mathbf {T}}^{\tfrac{1}{2}}_t \end{aligned}$$
(5)

where \(\exp _{\mathbf {X}}(\mathbf \Sigma )\) is an exponential mapping that maps a point \(\mathbf {X}\) on the manifold to another point on the manifold with a tangent vector \(\mathbf \Sigma \) and its inverse mapping can also be uniquely defined as a logarithmic mapping. A simple illustration is shown in Fig. 2. See [9] for more details.

The variance \(\sigma ^2\) is defined as

$$\begin{aligned} \sigma ^2 = \mathcal {E}[d(\bar{\mathbf {T}}, \mathbf {T})^2] = \frac{1}{N}\sum ^N_{i=1}d(\bar{\mathbf {T}}, \mathbf {T}_i)^2 \end{aligned}$$
(6)

where \(\mathcal {E}[\cdot ]\) denotes the expectation of a random tensor \(\mathbf {T}\). The distance, mean value and variance in a Riemannian metric space are utilized in a GMM to estimate the probability distribution of each component described in Sect. 2.3.

Fig. 1.
figure 1

A part of tensor field of CT slice. Tensors only with \(\mathcal {V} > 0\) are illustrated.

Full size image
Fig. 2.
figure 2

Illustration of manifold mapping. \(\mathbf {x}\) and \(\mathbf {y}\) are two tensors (matrices) on the manifold \(\mathcal {M}\) [11]. \(\varvec{\Sigma }\) is a tensor on the tangent space \(T_{\mathbf {x}}\mathcal {M}\) at tensor \(\mathbf {x}\). The exponential mapping \(\exp (\cdot )\) projects each point at the tangent space onto the manifold

Full size image

2.3 Graph-Cut Utilizing a Tensor Riemannian Metric

Graph-cut is a powerful optimization algorithm proposed by Boykov et al. [12] in 2001. The minimal cut corresponding to a optimization solution of an energy function E can be found effectively based on the max-flow/min-cut theorem. Whereas the conventional graph-cut segmentation method only utilizes low-order information such as intensity, the proposed tensor-based graph-cut method makes use of high-order information, the local geometric structure. The energy function in this work is given by

$$\begin{aligned} \begin{aligned} E(L)&= \underbrace{\underbrace{\omega \sum _{\mathbf {x}\in \mathbb {X}}-\log \Pr (\mathbf {x}|L_{\mathbf {x}})}_{\text {intensity term}} + \underbrace{(1 - \omega )\sum _{\mathbf {T}\in \mathbb {T}}-\log \Pr (\mathbf {T}|L_{\mathbf {T}})}_{\text {tensor term}}}_{\text {data term}}\\&\quad + \underbrace{ \varphi {\sum _{\{\mathbf {x}_m, \mathbf {x}_n\}\in \mathcal {N}}}V_{m,n}(\mathbf {x}_m, \mathbf {x}_n) + (1 - \varphi ){\sum _{\{\mathbf {T}_m, \mathbf {T}_n\}\in \mathcal {N}'}}U_{m,n}(\mathbf {T}_m, \mathbf {T}_n)}_{\text {smoothness term}}\\ \end{aligned} \end{aligned}$$
(7)

where \(\omega \) and \(\varphi \) are the weight parameters between the intensity-based term and the tensor-based term. \(\Pr (\cdot )\) is the GMM probability that the voxel \(\mathbf {x} \in \mathbb {X}\) or tensor \(\mathbf {T} \in \mathbb {T}\) be assigned to the label \(\{L_\text {B}, L_\text {F}\}\), \(L_\text {B}\) and \(L_\text {F}\) denote the background label and foreground label, respectively. Here, \(\mathbf {x}_m\) and \(\mathbf {x}_n\), as well as \(\mathbf {T}_m\) and \(\mathbf {T}_n\) are two pairs of neighboring voxels belonging to neighboring pair sets \(\mathcal {N}\) and \(\mathcal {N}'\). Finally, \(V(\cdot , \cdot )\) and \(U(\cdot , \cdot )\) represent the dissimilarity measures for neighboring voxels and neighboring tensors, respectively. Here, we focus on the tensor term. The GMM distribution can be defined as \(\Pr (\mathbf {T}|\bar{\mathbf {T}}, \sigma ) = \sum _k\eta _k\Pr (\mathbf {T}|\bar{\mathbf {T}}_k,\sigma _k)\), where k denotes the k-th component of GMM, and \(\eta _k\) is the mixture weight parameter. Therefore, the k-th Gaussian distribution of tensors can be formulated as:

$$\begin{aligned} \Pr (\mathbf {T}|\bar{\mathbf {T}}_k, \sigma _k) = \frac{1}{\sigma _k\sqrt{2\pi }}\exp \left( -\frac{d^2(\mathbf {T}, \bar{\mathbf {T}}_k)}{2\sigma _k^2}\right) . \end{aligned}$$
(8)

The dissimilarity measure of tensors is given by

$$\begin{aligned} U_{m,n}(\mathbf {T}_m, \mathbf {T}_n) = \frac{\exp \left( -\xi d^2(\mathbf {T}_m \mathbf {T}_n)\right) }{dist(\mathbf {T}_m, \mathbf {T}_n)} \end{aligned}$$
(9)

where \(\xi \) is a contrast adaptive constant, and \(dist(\cdot , \cdot )\) is the Euclidean distance.

In this work, the initial user-specified “trimap” is generated automatically. K-means clustering algorithm is applied to the vesselness enhancement result \(I_{\mathcal {V}}\) to extract the most probable blood vessel regions as foreground regions. The procedure can then be written as:

$$\begin{aligned} \mathbf {C}'&= \mathop {\hbox {argmin}}\limits _{\mathbf {C}}\sum _{i=1}^{K}\sum _{\mathbf {x}\in \mathbf {c}_i}||I_{\mathcal V}(\mathbf {x}) - \bar{I}_i ||^2 \end{aligned}$$
(10)
$$\begin{aligned} \mathbf {c}'_{max}&= \mathop {\hbox {argmax}}\limits _{\mathbf {c}'_i}\frac{1}{|\mathbf {c}'_i|}\sum _{\mathbf {x}\in \mathbf {c}'_i}|I_{\mathcal V}(\mathbf {x})| \end{aligned}$$
(11)

where the set \(\mathbf {C}\) consists of K initial clusters \(\mathbf {c}_i\) \((i=1,\dots ,K)\) and \(\mathbf {C}'\) is the final clusters. \(\bar{I}_i\) is the average intensity of set \(\mathbf {c}_i\), \(|\mathbf {c}'_i|\) is the number of the voxels in i-th cluster \(\mathbf {c}'_i\). Voxel set \(\mathbf {c}'_{max}\) is considered as most probable blood vessel regions. This scheme can effectively reduce the user’s interaction complexity.

3 Experiments and Results

The proposed method was tested on seven cases of contrast-enhanced renal CT data, including six cases of renal cancer and one normal case, taken in arterial phase. All the experiments were performed on a VOI of kidney, which is a cubical data created manually containing the whole kidney region. The size, pixel spacing as well as slice spacing of the VOI data were set to \(150\times 150\times 200\) voxels, 0.65–0.70 mm and 0.40–0.50 mm, respectively. The gold standards were created by two human experts having medical knowledge. We implemented the algorithm in C++ and evaluated it on a normal computer with an Intel Xeon CPU E5-2667 2.9 GHz 6 cores. Average computation time for one case is about 10 min. In Eq. 7, the weight parameters \(\omega \) and \(\phi \) are set to 0.8–0.9 through all experiments. Minimum and maximum detect size of Hessian vesselness filter are 1 mm and 2 mm, detect step is 1 mm.

Fig. 3.
figure 3

Experimental results of three cases. Each row represents one single case. (a) is gold-standard data, (b) is the result of graph-cut in previous work [5], (c) is the result of proposed graph-cut method. (d) is a part of CT images. The red circles correspond to the renal arteries marked in the (c). (Color figure online)

Full size image

In this study, four validation metrics, i.e. Dice Coefficient (DC), Overlapping of Centerline (OC), True Positive Ratio (TPR) and False Positive Ratio (FPR), were utilized to evaluate the proposed algorithm. OC index represents segmentation accuracy in the meaning of the branch length and the number of branches. Also, DC is the metric based on regions that is easily affected by thick branches. Both experimental results for different cases and the comparison between the previous and the new proposed graph-cut methods have been shown in Fig. 3. Quantitative evaluations of the proposed method versus the previous method have been presented in Table 1.

Table 1. Quantitative results of proposed method. For each case, boldface numerals show results of proposed method and lightface numerals show results of previous method [5]
Full size table

4 Discussion and Conclusions

We have presented a tensor-based graph-cut method for blood vessel segmentation in renal CT. Firstly, a tensor field is generated by calculating the Hessian matrix. Secondly, Riemannian metric is utilized to compute statistics over tensors. Finally, these statistics are used in graph-cut framework to estimate the probability distribution of foreground and background regions. According to the OC index, it is evident that the tensor-based graph-cut has an excellent segmentation performance at tiny blood vessels. It can extract more tiny blood vessels than previous method that is a hybrid method involving graph-cut and vessel tracking techniques. In Case 5, the previous method has a higher OC index than the proposed method, but also has a higher FPR. Serious over-segmentation is also shown in Fig. 3.

We used both intensity and tensors in the energy function of Eq. (7). This is because “tensors” only contain geometric information. However, we can only distinguish arteries from veins by the intensity difference. Thus, we solved this problem through a linear combination of intensity term and tensor term. But high weight values for tensor term may lead to over-segmentation by growing the segmentation into other tubular structures. This is the reason that the proposed method has a lower DC than previous method in several cases. For future work, we are going to find an automated method to decide the optimal weight parameters \(\omega , \varphi \). In addition, more experiments will be tested to assess the robustness of the proposed algorithm.