Automatic coronary isolation using a n-MIP ray casting technique

[0001] This application claims priority to U.S. Provisional Application No. 60/520,789, which was filed on Nov. 17, 2003, and which is fully incorporated herein by reference.

[0002] 1. Field of the Invention

[0003] The present invention relates generally to the field of medical imaging, and, more particularly, to automatic coronary isolation using a n-MIP ray casting technique.

[0004] 2. Description of the Related Art

[0005] Heart disease is a leading cause of death in the United States. One type of heart disease is coronary artery disease. Coronary artery disease is a narrowing or blockage of coronary arteries, which are the arteries that provide the heart with blood. This narrowing or blockage potentially slows blood flow, which decreases the oxygen and nutrients going to the heart. This may lead to chest pain (i.e., angina), or, if one or more of the coronary arteries are completely blocked, a heart attack. Other problems resulting from coronary pathologies include congestive heard failure (hereinafter “CHF”).

[0006] Coronary arteries become narrowed or clogged by plaque build-up in the artery walls. Plaque may include excessive cholesterol and other substances in the bloodstream, such as inflammatory cells, proteins and calcium. Plaque build-up in the arteries is called atherosclerosis, which is also known as “hardening of the arteries.”

[0007] The early detection of plaque, aneurysms and abnormal configuration of coronary arteries is an important factor in the diagnosis and treatment of coronary pathologies. Echocardiography is a leading tool for the detection of plaque, aneurysms and abnormal configuration of coronary arteries. An echocardiogram uses high-frequency sound waves to visualize the anatomy, structure, and function of the heart. Echocardiography disadvantages include low resolution and noise. Generally, it is not possible to visualize coronary arteries in echocardiography.

[0008] More recently, computer tomography (“CT”) has increased in popularity in a wide variety of medical diagnosis fields, such as the detection of cancer. CT provides sub-millimetric spatial resolution and low noise level. However, a very limited number of techniques exist for the detection of heart peripheral vessels in CT.

[0009] Magnetic Resonance (“MR”) provides the same resolution and low noise level provided by CT, and, additionally, is totally non-invasive. One using MR does not need any contrast agent injection as required by CT. Furthermore, CT uses X-rays, which are potential health risks over long time exposure, while MR uses the body's magnetic field. However, users of pacemakers, aneurysm clips, and other implants cannot use MR.

[0010] In one aspect of the present invention, a method of automatically retrieving centerlines of peripheral vessels of a heart in an image is provided. The method includes segmenting the heart in the image to produce a three-dimensional (“3D”) segmented myocardium; computing a 3D distance map from the three-dimensional segmented myocardium; casting normal maximum intensity projection rays for each voxel on the surface of the heart in the 3D distance map; filtering the rays to obtain vessel candidate points; and refining the vessel candidates points to obtain centerlines.

[0011] In another aspect of the present invention, a method of automatically retrieving centerlines of peripheral vessels of a heart in an image is provided. The method includes segmenting the heart in the image to produce a three-dimensional (“3D”) segmented myocardium; computing a 3D distance map from the three-dimensional segmented myocardium; modeling a heart wall by a spheroid using the 3D distance map; casting normal maximum intensity projection rays for each voxel on the spheroid; filtering the rays to obtain vessel candidate points; and refining the vessel candidates points to obtain centerlines.

[0012] In yet another aspect of the present invention, a machine-readable medium having instructions stored thereon for execution by a processor to perform method of automatically retrieving centerlines of peripheral vessels of a heart in an image is provided. The method includes segmenting the heart in the image to produce a three-dimensional (“3D”) segmented myocardium; computing a 3D distance map from the three-dimensional segmented myocardium; casting normal maximum intensity projection rays for each voxel on the surface of the heart in the 3D distance map; filtering the rays to obtain vessel candidate points; and refining the vessel candidates points to obtain centerlines.

[0013] In a further aspect of the present invention, a machine-readable medium having instructions stored thereon for execution by a processor to perform method of automatically retrieving centerlines of peripheral vessels of a heart in an image is provided. The method includes segmenting the heart in the image to produce a three-dimensional (“3D”) segmented myocardium; computing a 3D distance map from the three-dimensional segmented myocardium; modeling a heart wall by a spheroid using the 3D distance map; casting normal maximum intensity projection rays for each voxel on the spheroid; filtering the rays to obtain vessel candidate points; and refining the vessel candidates points to obtain centerlines.

[0021] Illustrative embodiments of the invention are described below. In the interest of clarity, not all features of an actual implementation are described in this specification. It will of course be appreciated that in the development of any such actual embodiment, numerous implementation-specific decisions must be made to achieve the developers' specific goals, such as compliance with system-related and business-related constraints, which will vary from one implementation to another. Moreover, it will be appreciated that such a development effort might be complex and time-consuming, but would nevertheless be a routine undertaking for those of ordinary skill in the art having the benefit of this disclosure.

[0022] While the invention is susceptible to various modifications and alternative forms, specific embodiments thereof have been shown by way of example in the drawings and are herein described in detail. It should be understood, however, that the description herein of specific embodiments is not intended to limit the invention to the particular forms disclosed, but on the contrary, the intention is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the invention as defined by the appended claims.

[0023] It is to be understood that the systems and methods described herein may be implemented in various forms of hardware, software, firmware, special purpose processors, or a combination thereof. In particular, at least a portion of the present invention is preferably implemented as an application comprising program instructions that are tangibly embodied on one or more program storage devices (e.g., hard disk, magnetic floppy disk, RAM, ROM, CD ROM, etc.) and executable by any device or machine comprising suitable architecture, such as a general purpose digital computer having a processor, memory, and input/output interfaces. It is to be further understood that, because some of the constituent system components and process steps depicted in the accompanying Figures are preferably implemented in software, the connections between system modules (or the logic flow of method steps) may differ depending upon the manner in which the present invention is programmed. Given the teachings herein, one of ordinary skill in the related art will be able to contemplate these and similar implementations of the present invention.

[0024] We present a heart peripheral vessels reconstruction method that assumes the existence of a segmented volume representing the heart myocardium. The method is based upon the fact that the vessels are generally parallel to the heart surface. Therefore, segmenting the heart wall may provide important information about the blood vessels. The methods described herein can be used for MR and CT modality images.

[0025] We consider a novel, multi-stage approach for automatic coronary vessel segmentation.

[0026] (1) The heart is segmented, and the heart wall shell is acquired.

[0027] (2) The surface of the heart is modeled by a simple geometrical volume, such as a spheroid.

[0028] (3) A ray is cast from each point on the spheroid surface, forming a three-dimensional (“3D”) ray collection, and the intersected intensities are recorded. The 3D ray collection is used as a pre-segmentation tool. Each vessel crossed by a ray generates a peak of intensity on the profile curve of the ray. This is a simple technique to detect voxels belonging to vessels. High-intensity tubular structures in this voxel space can then be used to detect the vessels.

[0029] (4) A full vessel tree is built, using vessel tracking techniques and minimum spanning tree.

[0030] Heart Segmentation and Distance May Computation

[0031] We consider a segmentation method driven from a graph optimization technique with a shape constraint. The idea lying behind this graphcut segmentation is to minimize an energy function that is defined on a graph, according to the cut of minimum weight. The energy is written as the sum of two terms: E_smooth(f) that imposes smoothness constraints on the segmentation map, and E_data(f) measuring how the label f is adapted to the data:

[0032] $\begin{array}{cc}E\left(f\right)=E_{\mathrm{smooth}}\left(f\right)+E_{\mathrm{data}}\left(f\right),& \left(1\right)\\ E_{\mathrm{smooth}}\left(f\right)=\sum _{p,q\in \mathrm{neighbors}}^{}{V}_{p,q}\left(f\left(p\right),f\left(q\right)\right),& \left(2\right)\\ E_{\mathrm{data}}\left(f\right)=\sum _{p\in P}^{}{D}_p\left(f\left(p\right)\right).& \left(3\right)\end{array}$

[0033] V_p,q in equation (2) is the interaction function between the pair of neighboring pixels {p, q}. D_p in equation (3) measures how close the label f is to the pixel p intensity. Generally, such a method provides a global optimal solution for the case of binary valued f(p).

[0034] The heart segmentation method provided above is only exemplary. It should be appreciated that other methods for isolating the art may be used, as contemplated by those skilled in the art, such as a model-based segmentation and segmentation methods based on level set methods.

[0035] The segmentation described above produces a 3D mask (i.e., pixels labeled object and background). A distance map from the surface of the 3D mask can provide valuable constraints during the recovery of the peripheral vessels. For example, as the peripheral vessels stay parallel to the surface of the heart, their distance (i.e., distance of the vessels from the heart's surface) in the distance map varies smoothly.

[0036] The distance map is computed by parsing the 3D mask twice—once in one direction, and once in the other direction. A distance map is basically a volume that associates to each point the distance between this point and a reference (in our case, the heart's surface). Each voxel in the distance map is filtered on an edge (object-background) by a 3D chamfer mask, M, as is commonly known to those skilled in the art. The filtered voxels are used to model the heart by a simpler geometrical object, such as an ellipsoid, to flatten its surface easily using cartography methods, as described in greater detail below.

[0037] Modelization of the Shell by a Spheroid

[0038] The distance map is used to model the heart wall by an ellipsoid or a biaxial spheroid. Although one may consider a more accurate model such as a tri-axial ellipsoid, a major drawback is that no simple mathematical solution to the 3D-2D projection exists.

[0039] The biaxial ellipsoid projection is a well-known technique in cartography, but not in medical imaging of the heart The biaxial spheroid reduces the deformations potentially occurring if using a sphere. For a biaxial ellipsoid of semi-axes length a and b, the surface equation is

[0040] $\begin{array}{cc}\frac{x^2}{a^2}+\frac{y^2}{a^2}+\frac{x^2}{b^2}=1& \left(4\right)\end{array}$
or, in a parametric form: x=a cos(λ)sin(φ), y=a sin(λ)sin(φ), z=b cos(φ) where λε[0,2π] and φε[0,π]. Similar to what is commonly used in cartography, λ and φ are called longitude and latitude respectively. The ellipsoid center is computed as the center of gravity. G of all the points located on the distance map isosurface 0 (the heart shell). The large axis Z is the vector {right arrow over (GM)}, where M is the point on the isosurface maximizing the length ∥{right arrow over (GM)}∥. Similarly, the small axis X is the vector {right arrow over (GN)}, where N is the point on the isosurface minimizing the length ∥{right arrow over (GN)}∥. The axis Y is deduced to have a direct orthogonal base, B=(G, X, Y, Z).

[0041] It should be noted that the quality of the modeling does not rely on the quality of the segmentation, which makes the method described herein independent from the selection of the segmentation method described above, and robust to noise. Moreover, unlike other organs, the heart naturally has the shape of an ellipsoid. From the ellipsoid surface, rays are cast to compute a two-dimensional (“2D”) view of the heart surface.

[0042] Ray Casting and N-MIP Projection

[0043] Next, we cast rays from the ellipsoid, and collect the voxel intensities in a predefined direction and range inside and outside the heart wall. The distribution of the nodes on the ellipsoid used to cast rays is computed through the following simple transformation. Referring now to FIG. 1, for each point P(λ, ) 105 in a 3D-space (on an ellipsoid), the 2D-point X(u, v) 110 (on a plane) is computed according to

[0044] $\begin{array}{c}\left[-\pi ,\pi \right]\left[-\pi /2,\pi /2\right]\mapsto \left[0,w\right]\left[0,H\right],\\ \left(\lambda ,\varphi \right)\mapsto \left(u=\frac{W}{2\pi }\lambda ,v=\frac{H}{\pi }\varphi \right).\end{array}$

[0045] If rays are cast directly from the distance map itself, the rays would not be homogeneous enough to be exploitable. Indeed, because of the heart surface irregularities, such rays would miss many vessels. Another solution may be based on transverse Mercator projection for a biaxial ellipsoid. The drawback of such a method is that the volume has to be an ellipsoid of revolution, meaning the semi-minor axis is on the polar plane, and the semi-major axes are on the equatorial plane.

[0046] Referring now to FIG. 2, examples of unfolded n-MIP (i.e., normal maximum intensity projection) views of heart peripheral vessels for various patients are shown with inverted colors. As the intensities have been inverted for the views to look like angiograms, the coronary arteries appear darker as compared to the other heart components. Each pixel I(x, y) on the 2D view is the normalized sum of each intensity, I_n, of the ray

[0047] $\begin{array}{cc}I\left(x,y\right)=\frac{1}{N}\sum _{0<n<N+1}^{}\mathrm{In}& \left(5\right)\end{array}$

[0048] A profile matching technique is used to detect whether a ray crosses a vessel structure. Referring now to FIG. 3, the ray profile curve starts in the heart wall, which is represented as a large, flat, high intensity area. FIG. 3 shows profile curves at two different locations on the surface of the heart. The myocardium area is indicated by crosshatching. The vascular structures appear as peaks of high intensity, of limited diameter, and at a certain distance from the pericardium.

[0049] Next, peaks are processed to determine if they are suitable candidates for a vessel point, according to the following criteria

[0050] 1. It respects certain intensity properties (intensity value and peak shape)

[0051] 2. It is within a certain distance from the heart wall

The intensity peaks, added to the local maximum (a ridge detection), allow a fair detection of the vessels on the rays profile curve. To detect peaks, we use a zero-crossing of the Laplacian, with the following kernel: [−1 2 −1].

[0052] Referring now to FIG. 4, upon completion of such a procedure, the peaks detected on the rays are associated to 3D points. As shown in FIG. 4, these 3D points may not be exactly at the center of the vessels, but they are next to or within the vessel lumen. Such a condition has its origin in the ray-casting effect, as the rays are homogeneously distributed through space. The next step will center these candidate points, filter the noise out, and track the detected vessels.

[0053] Refining the Vessel Candidate Points to Find the Centerlines

[0054] Vessel candidate voxels can be assumed to be next to or within a vessels lumen. Nevertheless, as shown in FIG. 4, the voxels are not yet on the vessels centerline. Furthermore, also as shown in FIG. 4, the peak detection may be sensitive to noise. Therefore, additional processing is required to center these points and to eliminate the noise. A refinement operation includes an eigenvalues analysis for each 3D-point within its corresponding local intensity space neighbors. With

[0055] $\begin{array}{cc}{A}_{\mathrm{ij}}=\frac{{\partial }^{2}I}{\partial i\partial j},& \left(6\right)\end{array}$
where i and j are equal to x, y or z, and the image intensity function I. At point P in space, the Hessian matrix M is defined as
M(P)=[A_ij]_{[i=x . . . z, j=x . . . z]}.
As the matrix M is defined, symmetric and positive, the computation of the eigenvalues is straightforward. At the center of the vessels, the three eigenvalues λ₁, λ₂ and λ₃ verify the following equations

[0056] There are two principal directions given by the eigenvectors. The vector v1 associated with λ1 corresponds to the local vessel orientation, whereas V2 and v3 define the vessel tangential plane.

[0057] For each candidate point P, the intensity function on the tangential cut plane defined by {right arrow over (v₂)} and {right arrow over (v₃)} is filtered by Gaussian functions, G_σ, with increasing variance σ. The G_σ minimizing the L² norm of (I−G_σ) is kept as model. Then, a potential function V_σ is built and minimized. V_σ represents how well a Gaussian vessel model fits for a variance_σ at point P

[0058] $\begin{array}{c}A=\frac{{\lambda }_2}{{\lambda }_3}\\ B=\frac{{\lambda }_2}{\sqrt{{\lambda }_3{\lambda }_2}}\\ S=\sqrt{{\lambda }_1+{\lambda }_2+{\lambda }_3}\\ V_{\sigma }\left(p\right)=\left(1-e^{\frac{A^2}{2{\alpha }^2}}\right)e^{\frac{B^2}{2{\beta }^2}}\left(1-e^{\frac{A^S}{2{\gamma }^2}}\right)+\int {G_{\sigma }-I}^2\end{array}$
where α, β and γ are used to control the sensitivity of V_σ to the diverse ratio A, B and S. This approach is based on an intensity analysis. Other approaches would include, for instance, taking the mean of the candidate points over a neighborhood. However, noise and high spacing between the points depending on the density of rays make these methods unsuitable for the present application. To account for the high spacing between vessel points, we consider a vessel tracking technique to increase the number of points and link the points together to produce a complete coronary arteries tree.

[0059] After the eigenanalysis, we assume we have the points, correctly placed, even though their number is not sufficient to track a vessels centerline. Thus, to place more points, a simple vessel tracking operation is performed for each 3D candidate point, under the conditions of linear intensity variation, consistent distance from the segmented heart wall, and consistent orientation. In other words, from a point resulting from the eigenanalysis, we build another point and link it to the previous one, thus tracking the vessel. The speed vector {right arrow over (v_t)} is the weighted sum of two vectors. The speed vector {right arrow over (v_t−1)} and the local orientation vector {right arrow over (v₁)} from the local vessel orientation.
{right arrow over (v_t)}=λ{right arrow over (v_t−1)}+(1−λ){right arrow over (v₁)}
where λε[0,1]. This speed vector, as well as intensity and distance of the current candidate point to the heart, has to be homogeneous. We detected edges through a Laplacian analysis:
∫|{right arrow over (∇)}{right arrow over (dl)}|<_coronary−I_exterior.
If such an inequality is not satisfied for the distance function D and the image intensity I, an edge has been detected, and the tracking operation stops. Otherwise, this tracking operation is performed until we reach the end of the vessel, or until we can link all the chains to another point. At this stage, the candidate points are centered and homogeneously distributed on the centerline of the vessel. They are linked by a minimum spanning tree algorithm, minimizing the following cost function E between two points P₁ and P₂
E_P1P2=|D(P₁)−D(P₂)|²+|I(P₁)−I(P₂)|²+∥{right arrow over (P₁P₂)}∥²+|{right arrow over (v₁)}(P₁)□{right arrow over (v₁)}(P₂)|,
and result in the peripheral vessels centerline, as illustrated in FIG. 5. More specifically, FIG. 5 shows points on the vessels centerline linked together.

[0060] Referring now to FIG. 6, a heart peripheral vessel segmentation method based on a ray filtering method is shown, as described in greater detail above, and in accordance one exemplary embodiment of the present invention. In a step-by-step approach, the heart is first segmented (at 605) using an method based on graph cuts and geodesic surfaces, for example. A 3D distance map is computed (at 610) out of the segmentations output. Then, from this distance map, the biaxial ellipsoid modeling the heart is computed (at 615). Rays are cast (at 620) from the surface of the ellipsoid toward the heart. Then, candidate points likely belonging to vessels are detected (at 625) along the n-MIP rays profiles. Finally, the vessels centerline is tracked (at 630) from the candidate points.

[0061] This method can be used, for example, to visualize the full peripheral vessel tree, and detect plaques and aneurysms. Once segmented, the vascular structures can be unfolded, and the quantification of the stenosis and the aneurysms is straightforward. The segmentation results support efficient reporting by enabling automatic generation of overview visualizations, guidance for virtual endoscopy, generation of curved MPRs along the vessels, or cross-sectional area graphs. Moreover, as applied to CT data sets, this method detects vascular structures quickly enough for industrial applications (60 seconds for a bi-processor 900 MHz, 1 GB RAM). Further, this method, as compared to techniques relying on front propagation, is robust to noise.

[0062] The particular embodiments disclosed above are illustrative only, as the invention may be modified and practiced in different but equivalent manners apparent to those skilled in the art having the benefit of the teachings herein. Furthermore, no limitations are intended to the details of construction or design herein shown, other than as described in the claims below. It is therefore evident that the particular embodiments disclosed above may be altered or modified and all such variations are considered within the scope and spirit of the invention. Accordingly, the protection sought herein is as set forth in the claims below.