3-D reconstruction of the coronary artery tree from multiple views of a rotational X-ray angiography

Abstract

To present an efficient and robust method for 3-D reconstruction of the coronary artery tree from multiple ECG-gated views of an X-ray angiography. 2-D coronary artery centerlines are extracted automatically from X-ray projection images using an enhanced multi-scale analysis. For the difficult data with low vessel contrast, a semi-automatic tool based on fast marching method is implemented to allow manual correction of automatically-extracted 2-D centerlines. First, we formulate the 3-D symbolic reconstruction of coronary arteries from multiple views as an energy minimization problem incorporating a soft epipolar line constraint and a smoothness term evaluated in 3-D. The proposed formulation results in the robustness of the reconstruction to the imperfectness in 2-D centerline extraction, as well as the reconstructed coronary artery tree being inherently smooth in 3-D. We further propose to solve the energy minimization problem using α-expansion moves of Graph Cuts, a powerful optimization technique that yields a local minimum in a strong sense at a relatively low computational complexity. We show experimental results on a synthetic coronary phantom, a porcine data set and 11 patient data sets. For the coronary phantom, results obtained using different number of views are presented. 3-D reconstruction error evaluated by the mean plus one standard deviation is below one millimeter when 4 or more views are used. For real data, reconstruction using 4 to 5 views and 256 depth labels averaged around 12 s on a computer with 2.13 GHz Intel Pentium M and achieves a mean 2-D back-projection error of 1.18 mm (ranging from 0.84 to 1.71 mm) in 12 cases. The accuracy for multi-view reconstruction of the coronary artery tree as reported from the phantom and patient studies is promising, and the efficiency is significantly improved compared to other approaches reported in the literature, which range from a few to tens of minutes. Visually good and smooth reconstruction is demonstrated.

Appendix

A. Coherence enhanced diffusion

For CED, The diffusion equation applied on the image I is given by:

$$ \partial_{t} I(x,t) = \nabla \cdot (D(\nabla I(x,t))), $$

(7)

where $\nabla \cdot$ is the divergence operator, $D$ is a diffusion tensor that depends on local image structures so that not only the amount, but also the direction of the diffusion can be regulated. Let $K_{\sigma }$ be the Gaussian kernel:

$$ K_{\sigma } (x) = {\frac{1}{{2\pi \sigma^{2} }}}e^{{ - {\frac{{\left| x \right|^{2} }}{{2\sigma^{2} }}}}} $$

(8)

and $I_{\sigma }$ be the convolving result of $I$ with $K_{\sigma }$: $I_{\sigma } (x,t) = K_{\sigma } * I(x,t)$.

The structure tensor $J_{\rho }$ is defined as the convolving result of gradient tensor product with $K_{\rho }$:

$$ J_{\rho } (\nabla I_{\sigma } ) = K_{\rho } * (\nabla I_{\sigma } \otimes \nabla I_{\sigma } ) . $$

(9)

$J_{\rho }$ has the eigenvalues $\mu_{1}$ and $\mu_{2}$, and the corresponding eigenvectors $w_{1}$ and $w_{2}$. The coherence-enhancing anisotropic diffusion tensor D is structured to have the same eigenvectors as $J_{\rho }$, but with customized eigenvalues $\lambda_{1}$ and $\lambda_{2}$,

$$ D = \left( {w_{1} |w_{2} } \right)\left( {\begin{array}{*{20}c} {\lambda_{1} } & 0 \\ 0 & {\lambda_{2} } \\ \end{array} } \right)\left( {\begin{array}{*{20}c} {w_{1}^{T} } \\ {w_{2}^{T} } \\ \end{array} } \right) $$

(10)

where

$$ \begin{aligned} \lambda_{1} &=& \kappa \\ \lambda_{2} &=& \kappa + (1 - \kappa )e^{{{\frac{ - 1}{{(\mu_{1} - \mu_{2} )^{2} }}}}} \\ \end{aligned} $$

(11)

with $\kappa \in (0,1)$. The diffusion process acts preferably along $w_{2}$ direction of the highest coherence within the scale $\rho$. $\sigma = 1.0$, $\rho = 4.0$, $\kappa = 0.001$, $t = 3$ are adopted in our algorithm.

B. Multi-scale hessian-based vesselness filtering

The MSHBVF is defined as:

$$ F(u) = \max\nolimits_{\varepsilon } f(u,\varepsilon ) $$

(12)

where f is the filter at scale ɛ. The filter response is calculated for different scales and the maximal response is used as the final vesselness result. The filter calculates 2^nd-order derivatives to build the Hessian matrix:

$$ \user2{H}_{\varepsilon } = \varepsilon^{2} \left[ \begin{gathered} \begin{array}{*{20}c} {{\frac{{\partial^{2} {\mathbf{u}}_{\varepsilon } }}{{\partial x^{2} }}}} & {{\frac{{\partial^{2} {\mathbf{u}}_{\varepsilon } }}{\partial x\,\partial y}}} \\ \end{array} \hfill \\ \begin{array}{*{20}c} {{\frac{{\partial^{2} {\mathbf{u}}_{\varepsilon } }}{\partial x\,\partial y}}} & {{\frac{{\partial^{2} {\mathbf{u}}_{\varepsilon } }}{{\partial y^{2} }}}} \\ \end{array} \hfill \\ \end{gathered} \right]. $$

(13)

Here ${\mathbf{u}}_{\varepsilon }$ is the convolution of the CED preprocessed image with a Gaussian kernel of size ε corresponding to the radius of the vessel of interest. Based on the hypotheses that the dark tubular structures on the bright background have the eigenvalues $\nu_{1} \approx 0$, $\nu_{2} > 0$, and $|\nu_{1} | \ll |\nu_{2} |$, the vesselness strength is calculated as [28]:

$$ V_{\varepsilon } = \left\{ {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 \\ {e^{{{{ - R_{B}^{2} } \mathord{\left/ {\vphantom {{ - R_{B}^{2} } {2\beta^{2} }}} \right. \kern-\nulldelimiterspace} {2\beta^{2} }}}} (1 - e^{{{{ - S^{2} } \mathord{\left/ {\vphantom {{ - S^{2} } {2\gamma^{2} }}} \right. \kern-\nulldelimiterspace} {2\gamma^{2} }}}} )} \\ \end{array} } & {\begin{array}{*{20}c} {\nu_{2} \le 0} \\ {\nu_{2} > 0} \\ \end{array} } \\ \end{array} } \right\} $$

(14)

where $R_{B} = \nu_{1} /\nu_{2}$ to enhance the tubular structure, $S = \sqrt {\nu_{1}^{2} + \nu_{2}^{2} }$ to suppress the background noise, and the vesselness orientation is calculated as the eigenvector corresponding to $\nu_{1}$. Here β and γ control the sensitivity of the filter and are set to 1 and 3, respectively in our algorithm.

C. α-Expansion move of graph-cuts optimization

GCs optimization using α-expansion move is based on computing a labeling corresponding to a minimum cut on a graph $g_{\alpha } = \left\langle {V_{\alpha } ,E_{\alpha } } \right\rangle$. The structure of the graph is determined by the current partition and by the current label α. The set of vertices:

$$ V_{\alpha } = \left\{ {\alpha ,\bar{\alpha },P,\mathop \cup \limits_{\begin{subarray}{l} \{ p,q\} \in N \\ f_{p} \ne f_{q} \end{subarray} } a_{{\{ p,q\} }} } \right\} , $$

(15)

includes two terminals α and $\bar{\alpha }$, the centerline pixels P, and the auxiliary nodes $a_{{\{ p,q\} }}$ between neighboring centerline pixels p and q such that $f_{p} \ne f_{q}$. Each centerline pixel $p \in P$ is connected to the terminals α and $\bar{\alpha }$ by t-links $t_{p}^{\alpha }$ and $t_{p}^{{\bar{\alpha }}}$, respectively. Each pair of neighboring centerline pixels $\{ p,q\} \in N$ with $f_{p} = f_{q}$ is connected by an n-link $e_{{\{ p,q\} }}$. For each pair of neighboring centerline pixels $\{ p,q\} \in N$ with $f_{p} \ne f_{q}$, a triplet edges $E_{{\{ p,q\} }} = \{ e_{{\{ p,a\} }} ,e_{{\{ a,q\} }} ,t_{a}^{{\bar{\alpha }}} \}$ are created, where $a = a_{{\{ p,q\} }}$ is the auxiliary node, $e_{{\{ p,a\} }}$ and $e_{{\{ a,q\} }}$ connect pixels p and q and the t-link $t_{a}^{{\bar{\alpha }}}$ connects the auxiliary node $a_{{\{ p,q\} }}$ to the terminal $\bar{\alpha }$. So all the edges are:

$$ E_{\alpha } = \left\{ {\mathop \cup \limits_{p \in P} \left\{ {t_{p}^{\alpha } ,t_{p}^{{\bar{\alpha }}} } \right\},\mathop \cup \limits_{\begin{subarray}{l} \left\{ {p,q} \right\} \in N \\ f_{{p \ne f_{q} }} \end{subarray} } E_{{\{ p,q\} }} ,\mathop \cup \limits_{\begin{subarray}{l} \left\{ {p,q} \right\} \in N \\ f_{{p = f_{q} }} \end{subarray} } e_{{\{ p,q\} }} } \right\}. $$

(16)

Please refer to [21] for the weights assigned to the edges. After partition, the global solution of one α-expansion move is as follows: a centerline pixel p is assigned the label α if the cut seperates p from the terminal α, and p retains its old label if the cut separates p from the terminal $\bar{\alpha }$.