Abstract

I. Introduction

The objective is to build a 3-D dynamic model of the heart from patient data acquired in Multi-slice Spiral Computed Tomography (MSCT). This 3-D model is obtained from a dynamic volume sequence and includes cardiac cavities, arteries and veins extracted on each volume. The interest of this model is manifold: (1) evaluate the performances of the segmentation, reconstruction and registration methods involved in the automatic characterization of the pathologies, (2) assess the kinetic properties of the structures in order to evaluate the myocardium viability and perfusion, (3) offer a numerical atlas to the student or clinician for the anatomy learning, (4) provide a virtual surgery training tool.

In our case, we want to simulate the data acquisition from a rotational angiography system to test 3-D reconstruction algorithms of vascular structures. It consists of modeling a continuous rotation of the C-arm system around the 3-D dynamic model and to compute a series of projected images. A wide range of projections in the axial, caudal and cranial angulations can then obtained which will be used for the 3-D reconstruction of the coronary structures.

This paper presents the construction of this 3-D dynamic model. The acquisition of a dynamic volume sequence leads to very large data sets (several hundreds of megabytes) from which the structures must be extracted under very strict time constraints to face the clinical requirements. We already dealt with the cardiac cavities in a previous work [1]. The current objective is to extract the vascular structures in order to complete the 3-D model of the heart. Most of the current segmentation methods bring already relevant delineations of vessels but with too high computation time. They rely on intensity-based methods, generalized cylinder approximations, multiscale and skeletonization schemes, and deformable model approaches applied on successive 2-D slices or on volume data [2].

The challenge is therefore to reach similar performance in few minutes on a standard PC platform, by focusing the extraction process only inside the vascular structure. The present work relies on a geometric moment operator [3] and the response of a multiscale Hessian filter [4] to track a vessel, locally estimate the local diameter and orientation and detect the bifurcation.

The outline of this paper is as follows: section II describes the adopted methodology for the coronary arteries extraction. Section III presents some preliminary results and section IV concludes on the perspectives.

II. Methodology

The tracking algorithm relies on a local modeling of the vessel by a cylinder in a 3-D homogeneous space. The parameters of this cylinder (location of the center of gravity, and diameter) are estimated using a 3-D geometrical moment operator. Its main advantage is that it provides analytical expressions for the computation of these parameters. A multiscale filter based on eigenvalue analysis of the Hessian matrix is then applied to highlight tubular structures and coping with varying widths. It is also used to locally determine the principal direction of the vessel. A detection of possible bifurcation is performed then in the estimated direction. If one bifurcation is detected, some seed points are extracted which are saved in a list to be further taken as initial point of the tracking process. The next point P_i+1 on the centerline is searched for in the estimated direction at the current point P_i using the moment operator.

The tracking of a branch can be divided into six stages:

1. Interactive selection of a seed point P_{i = 0};

2. Position refinement of the point P_i

3. Estimation of the local diameter d_i vessel intensity I_v and background I_b;

4. Local direction estimation;

5. Bifurcation detection;

6. if a bifurcation is detected then search for some candidate points. These points will be stored in the seed point list for further tracking;

   else search for the next point P_i+1 in the estimated direction;

7. Go to (b) until stop criterion is satisfied.

The tracking process was performed for every seed point of the list, the latter one being updated during the tracking when a bifurcation was detected. A flowchart of the algorithm is given Fig. 1.

Open in a separate window

Fig. 1

Flowchart of the algorithm

2.1 Position refinement of P_i and vessel local diameter estimation

Starting from the hypothesis that the vessel can be locally modeled by a cylinder of center of gravity P_i, diameter d_i and orientation vector v_vessel in the 3-D space, we looked for estimating the center of gravity of this cylinder so that it corresponds to a point located on the central axis of the vessel. We applied the geometric central moments of order p+q+r to compute this local centroid and estimate the local diameter of the vessel. The coordinates of the local centroid (,, ) inside a spherical window is given by:

$$(\overline{x},\overline{y},\overline{z})=(\frac{M_{100}}{M_{000}},\frac{M_{010}}{M_{000}},\frac{M_{001}}{M_{000}})$$

(1)

The distance between the window center and the local centroid allows detecting whether the window is located on the center of a bright structure or not. We applied thus an iterative process to move the center of the window P′_i+1 towards the centroid P_i+1 till this centroid coincides with the center of this window (P′_i+1 = P′_i+1), in order to shift it to the central axis of the vessel (Fig. 2). We then estimated the local diameter of the vessel from the zero order moment M₀₀₀, the local intensity inside (I_v) and outside (I_b) the vessel:

Open in a separate window

Fig. 2

Iterative centering of the spherical window to move its center P′_i+1 towards its center of gravity Pi₊₁

$$d_i=2R{\left[1-{\left(\frac{3M_{000}/4\pi -I_v}{I_{b-}{I}_v}\right)}^{2/3}\right]}^{1/2}$$

(2)

where R is the radius of the spherical window. Optimal parameters were obtained when the size of the spherical window fitted the vessel size. The local diameter was thus estimated using a multi-resolution local moment computation [5].

The intensities I_v and I_b were updated at each step of the tracking to take into account the heterogeneity of the intravascular and environmental surroundings [5]. The background intensity I_b was set to the mean value computed outside the vessel in its normal direction to the currently estimated vessel axis and outside the vessel. In the same way, the vessel intensity I_v was equal to the mean intensity computed along the estimated vessel direction.

2.2 Local vessel direction estimation

A multiscale filter based on eigenvalue analysis of the Hessian matrix was locally applied at the estimated point P_i to enhance the tubular structure and estimate the orientation of the vessel. The response of this filter was computed at different scales σ and the eigenvalues were combined into a discriminant function to allow discerning tube-like, blob-like and plate-like structures [4–6].

Let H_σ(x) be the Hessian matrix defined at a given voxel x at scale σ:

$$H_{\sigma }(\mathbf{x})=\left[\begin{array}{ccc}{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{array}\right]$$

(3)

where I_αβ(x) denotes regularized derivatives of the image I(x), which are obtained by convolving the image using the Gaussian kernel G(x, σ) at scale σ,

$$I_{\alpha \beta }(\mathbf{x})={\sigma }^2\frac{{\partial }^{2}G(\mathbf{x},\sigma )}{\partial \alpha \partial \beta }\ast I(\mathbf{x})$$

(4)

$$G(\mathbf{x},\sigma )=\frac{1}{\sqrt{{(2\pi {\sigma }^2)}^3}}{e}^{\parallel \mathbf{x}\mid {\mid }^2/2{\sigma }^2}$$

(5)

Let λ₁, λ₂, λ₃ and v₁, v₂, v₃, be the eigenvalues and unit eigenvectors of H_σ(x).

The eigenvector v₁ in the direction of the vessel corresponds to the smallest eigenvalue λ₁, while the eigenvectors v₂ and v₃ form a base in the orthogonal plane to v₁ (Fig. 3). The relations between the eigenvalues are thus the following:

Open in a separate window

Fig. 3

Eigenvector approximate the local vessel direction

$$\begin{array}{l}\mid {\lambda }_1\mid \le \mid {\lambda }_2\mid \le \mid {\lambda }_3\mid \\ \mid {\lambda }_3\mid \approx \mid {\lambda }_2\mid >>\mid {\lambda }_1\mid \approx 0,\hspace{0.17em}\text{for\hspace{0.17em}a\hspace{0.17em}voxel\hspace{0.17em}belonging}\hspace{0.17em}\text{to\hspace{0.17em}the\hspace{0.17em}vessel}\end{array}$$

(6)

Frangi et al. [6] introduced two geometric ratios R_A, R_B in their vessel likeliness function V(x, σ):

$$V(\mathbf{x},\sigma )=\{\begin{array}{l}0,\hspace{0.17em}\text{if\hspace{0.17em}}{\lambda }_2>0\hspace{0.17em}\text{or\hspace{0.17em}}{\lambda }_3>0\\ \left(1-e^{-R_A^2/2{\alpha }^2}\right)e^{-R_B^2/2{\beta }^2}\left(1-e^{-S^2/2c^2}\right),\hspace{0.17em}\text{else}\end{array}$$

(7)

with

$$R_B=\frac{\mid {\lambda }_1\mid }{\mid {\lambda }_2{\lambda }_3\mid },R_A=\mid \frac{{\lambda }_2}{{\lambda }_3}\mid ,S=\sqrt{\sum _{j=1}^3{\lambda }_j}$$

(8)

The first ratio R_A was used to distinguish between a plate-like and a tube-like structure. The second one, R_B, addressed the deviation from a blob-like structure and S was a measure of second order structureness, which was used to reduce the response of the background voxels, α, β, c control the sensitivity of the filter to deviations in R_A, R_B and S. The response of the filter is expected to be maximum at a scale σ that approximates the radius of the vessel. The maximum value among the set of response computed at different scales was given by:

$$V(\mathbf{x})=\underset{{\sigma }_{min}\le \sigma \le {\sigma }_{max}}{max}V(\mathbf{x},\sigma )$$

(9)

The vessel direction v_vessel was given by the eigenvector v₁ at the scale σ that provided the optimal response of the filter.

III. Results

3.3 Experimental results

Coronary arteries were extracted on the 10 volumes of two sequences, focusing on the four main branches: the right coronary artery (RCA), the left anterior descending artery (LAD), the circumflex artery (CRX) and the first diagonal artery (DA). Fig. 4 illustrates the extraction results for one of the volume of one sequence. All the branches were extracted from only 2 interactively selected seed points. In this volume, most of the branches were correctly extracted. However, some problems may occur when the heart rate varies during the acquisition, inducing thus motion artifacts that lead to detection errors. A way to deal with that problem will be to use the information extracted on the other volumes and perform a motion estimation to recover the branch segment blurred by the motion in the current volume. Fig. 5 depicts the extracted network with estimated local diameter. Fig. 6 provides an estimation of the diameter along the LAD artery in the perpendicular planes to the vessel.

Open in a separate window

Fig. 4

Extraction results for one volume of the sequence. Only 2 seed points was interactively pointed. Slices are displayed in the axial (a), coronal (b) and sagittal (c) planes respectively. The extracted branches are superposed in color on the slices (violet: LAD, green: CX, red: DA and pink: RCA), d) Curvilinear MPR showing the CX artery

Open in a separate window

Fig. 5

Extracted 3D coronary tree from two points of view

Open in a separate window

Fig. 6

Diameter estimation for the CX artery (red dash line) and average diameter (blue solid line)

IV. Conclusion

An efficient model-based solution has been proposed for the 3-D tracking of vessels in MSCT volumes. It provides a first approximation of the vascular patterns and allows extracting the structure of the coronary network with a very few initial seed points (between 2 and 4 depending on the complexity of the structure) within 3 minutes on a PC P4-2.4G, 1G RAM. These preliminary results appear promising. Difficulties remain nevertheless in presence of motion artefacts. Further extensions will be to introduce statistical models to improve the coronary extraction as to exploit the other volumes of the sequence to deal with the motion artefacts.

References