WO2003071925A2

WO2003071925A2 - Refinement of isointensity surfaces

Info

Publication number: WO2003071925A2
Application number: PCT/US2003/005268
Authority: WO
Inventors: Peter J. Yim
Original assignee: The Government Of The United States Of America As Represented By The Secretary Of The Department Of Health And Human Services
Priority date: 2002-02-22
Filing date: 2003-02-21
Publication date: 2003-09-04
Also published as: AU2003215354A1; WO2003071925A3; AU2003215354A8; WO2003071925A8

Abstract

Beginning with an isosurface, image forces are determined for vertices of a mesh according to a three dimensional image of a tubular vessel, such as a blood vessel. The image forces are applied to relocate the vertices into relative positions corresponding to a contour of the vessel, thereby providing improved representations of vessels. The improved representations, for example, facilitate estimation of the degree of stenosis in blood vessels of subjects suffering from atherosclerotic disease.

Description

REFINEMENT OF ISOINTENSITY SURFACES

FIELD

This invention relates generally to mesh modeling techniques. More particularly, the invention relates to refining isosurfaces to model the luminal shape of blood vessels and for other medical purposes.

BACKGROUND

The simplest surface reconstruction obtained from an image is the isosurface. An isosurface is defined as a surface in which all the vertices are located at points in the image that share the same intensity.

f(x_v) = T Vv e S

Where x_v) is the linearly interpolated image intensty at the location of vertex v, Eis the iso-intensity value, and S is the set of vertices in the isosurface. The isosurface can, for example, be constructed by the Marching-Cubes algorithm [Lorensen, W.E. and Cline, H.E., "Marching Cubes: A High Resolution 3D Surface Construction Algorithm:, Computer Graphics (Proc.SIGGRAPH), 21(3): 163-16, 1987]. Isosurfaces of blood vessels, the colon, and other anatomical structures can have a highly realistic appearance. However, isosurfaces have limited value for shape analysis since the image intensity of an object and the background are usually inhomogeneous. Thus, a single iso-intensity value usually cannot provide an accurate segmentation of an entire anatomical structure. This effect is pronounced in magnetic resonance angiography (MRA). For example, spin dephasing at stenoses reduces the image intensity in the vicinity of stenoses [Mitsuzaki et al., "Delineation of simulated vascular stenosis with Gd-DTPA-enhanced 3D gradient echo MR angiography: an experimental study," J Comput Assist Tomogr, 24(1):77- 82, 2000]. As a result, isosurfaces of MRA tend to exaggerate the degree of stenosis. More sophisticated approaches to surface reconstruction use the gradient magnitude of the image to define object boundaries. Presumably object boundaries fall where the gradient magnitude is a peak or ridge. This principle has been incorporated in the Watershed algorithm, Level Set algorithm and in prior deformable models. The Watershed algorithm detects ridges in the image in a manner directly analogous to the definition of topological watersheds [C. Lantuejoul, "La squelettisation et son application aux mesures topologiques des mosaiques polycristallines," Theses Dr. Ing. Ecole des Mines de Paris, 1978; L. Vincent and P. Soille, "Watersheds in digital spaces: an efficient algorithm on immersion simulations," IEEE Transactions on Pattern Analysis and Machine

Intelligence, 13: 583-598, 1991]. The analogy arises from the conceptualization of the image as a relief map in which the elevations are determined by the image intensities. If the Watershed algorithm is applied to the gradient magnitude of the image, ridges correspond to boundaries of the object in the original image. This algorithm is most intuitive when applied to 2D images but extends directly to multidimensional images [Yim et al., "High-resolution four-dimensional surface reconstruction of the right ventricle and pulmonary arteries," SPIE Proceedings, 3038: 726-738, 1998]. The Watershed algorithm produces a segmentation of the image in the discrete domain which is not adequate for visualization. However, Watershed boundaries can be converted to surface meshes [Yim PJ, Summers RM, "Analytic surface reconstruction by local threshold estimation in the case of simple intensity contrasts," SPIE Proceedings, 3660: 288-300, 1999]. A limitation of the Watershed algorithm is that smoothing constraints cannot be imposed on the segmentation boundary. In the Level Set algorithm, a surface is also reconstructed by growth process

[Malladi, R., Sethian, J. A. and Vemuri, B. C," Shape Modeling with Front Propagation: A Level Set Approach," IEEE Transactions on Pattern Analysis and Machine Intelligence 17(2), 158-175. 1995]. A seed point or region is manually placed within an object and grows outwards at a rate that is inversely proportional to the gradient magnitude at every point in the boundary. The rate of growth of a segmentation boundary is also governed by the local curvature of the boundary. Constraints on the curvature are such that high-curvature protrusions are suppressed and the surface is effectively smoothed. The boundary of the growth region is implicitly defined as the level set of a function. The potential advantage of this methodology is that it requires minimal initialization. A disadvantage of this methodology is that no distinction can be made between real protrusions of the underlying object and protrusions or surface roughness due to image noise. Thus, constraints on the curvature of the boundary must be relaxed to reconstruct narrow sections of an object which entails a decrease in the smoothness of the surface. While this methodology has been shown to give adequate results when applied to vessels in 2D cross-sectional images of a normal vessel tree, further validity of this method has not been established [Wang et al., "Level Sets for Vascular Model

Construction in Computational Hemodynamics," Proceedings of ASME Summer Bioengineering Conference, Big Sky, Montana, pp. 721-722, 1999]. This method is similar in principle to deformable models with adaptive meshes [Mclnerney T, Terzopoulos D., "Topology adaptive deformable surfaces for medical image volume segmentation," IEEE Trans Med Imaging 18(10):840-50, 1999; Gill et al.,

"Accuracy and variability assessment of a semiautomatic technique for segmentation of the carotid arteries from three-dimensional ultrasound images," Medical Physics, 27:1333-1342, 2000].

A contrasting approach to that of Level Sets is that of deformable models with fixed topologies. In these methods, the surface mesh is initialized to reflect the known shape of the object. For vascular surface reconstruction, for example, a tubular configuration is appropriate since the vessel is typically tubular in shape. The initial location of the mesh is determined by semi-automated methods. Yim et al. [Yim et al., "Vessel Surface Reconstruction with a Tubular Deformable Model," IEEE Transactions on Medical Imaging, 20(12): 1411-1421, 2002] proposed a tubular deformable model for vascular surface reconstruction in which the axis of the vessel is defined manually. Based on the vessel axis a discrete tubular coordinate system is constructed with radial lines emanating from the axis at all discrete axial and circumferential locations. In this model, only the radial component of vertex position is smoothed in contrast to the more generic smoothing applied in the model of Frangi et al. [Frangi et al., "Model-based quantitation of 3-D magnetic resonance angiographic images," IEEE T. Med. Imaging, 18: 946-956, 1999]. By smoothing only the radial component of vertex position the model has no predisposition towards either contraction or dilation. The tubular deformable model incorporates a scheme for merging of radial lines to prevent intersections due to high curvature of the vessel axis. These models have been shown to produce reliable results for contrast-enhanced MRA of the carotid artery bifurcation. In that case, the deformable models are applied independently to each branch of the bifurcation and its continuation on the parent vessel. These deformable models have been found to be reliable for vascular regions such as the carotid bifurcation where the bifurcation is relatively symmetric; the branches have comparable sizes. However, tubular deformable models are less suitable for cases where the bifurcation is highly asymmetric. In those cases, such as the branching of the renal artery from the aorta, there is large discontinuity in the shape of the smaller branch vessel as it joins the parent vessel that cannot be accurately represented with a tubular model.

SUMMARY

An isosurface mesh is produced from a three dimensional image according to a selected intensity value. Vertices of the mesh are repeatedly relocated according to the composite of an image force, torsion force, and stretch force determined for each of the vertices, until a stopping point is reached. The image force for a vertex is determined according to the gradient of the gradient magnitude image evaluated at the vertex. The image force is applied to relocate the vertex into a relative position corresponding to a contour of the vessel.

An effective force on each vertex is produced by a torsion on triangles associated with the vertex. The torsion on each triangle rotates the triangle so as to decrease the difference between the normal direction of the triangle and the average of the normal direction of the neighboring triangles. In one embodiment, the torsion on a triangle is the cross-product of the normal vector to a triangle with the unit vector in the direction of the average of the normal vectors of the neighboring triangles. The effective force on a given vertex due to the torsion on one of the associated triangles is proportional to the cross-product between the moment arm from the vertex to the axis-of-rotation of the triangle and the torsion arm. The stretch force is an attractive force that acts to regularize the spacing between adjacent vertices to create an image having a uniform appearance, and is proportional to the distance between vertices. The attractive force between vertices, in some embodiments, is proportional to a component of the displacement between vertices that is tangent to the surface, thereby preventing shrinkage of the surface due to the attraction between vertices.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flow chart illustrating a method embodiment to process an isosurface.

FIG. 2 is an illustration of a neighborhood of triangles and vertices of a mesh embodiment.

FIG. 3 is an illustration of a triangle of a mesh embodiment.

FIG. 4 is an embodiment of a formula chart for computing forces on vertices of a mesh.

FIG. 5 is a block diagram of an apparatus to process an isosurface.

FIG. 6 is a series of images showing the surface reconstruction using the disclosed methods of a digital phantom of a hypothetical artery having a prescribed tubular geometry with a 70% stenosis. FIG. 7 is a series of anterior images showing a contrast-enhanced MRA image, surface reconstructions, and a corresponding X-ray angiographic image of the abdominal aorta and the left renal artery of a subject with mild stenosis at the origin of the left renal artery.

FIG. 8 is a series of anterior images showing a contrast-enhanced MRA image and surface reconstructions in a subject with bilateral atherosclerotic disease of the renal arteries.

FIG. 9 is a series of images showing a contrast-enhanced MRA image and surface reconstructions of a carotid artery with atherosclerotic disease. DESCRIPTION

In the following figures and description, like numbers refer to like elements. References to "one embodiment" or "an embodiment" do not necessarily refer to the same embodiment, although they may. In one aspect, an isosurface model of a luminal surface is subjected to vertex relocation. The vertex relocation may lead to a vascular model with more consistency and accuracy in identifying stenoses and other features the luminal surface.

FIG. 1 is a flow chart of a method embodiment 100 to process an isosurface. At 102 forces are determined to apply to the vertices of the surface. These forces may be applied to effect relocation of the vertex coordinates, in effect altering the contour of the surface. At 104 the vertices are relocated according to the forces. Relocation may be incremental, in proportion to the composite of the determined forces, and in a direction dictated by the forces. If at 106 the stopping point is reached, the method 100 concludes. Otherwise, forces are determined again at 102 for the new positions of the vertices and the method 100 repeats.

In one embodiment, the stopping point is reached after a preselected number of iterations of the method 100. In another embodiment, the stopping point is reached when the vertex relocations between successive iterations of the method 100 become insubstantial (e.g. successive iterations don't result in substantial changes in the positions of the vertices). The stopping point may be any selected threshold level below which relocation is regarded as insubstantial.

In one embodiment, each vertex of the mesh comprises a three dimensional coordinate. The forces on a particular vertex reflect the amount and direction of relocation to apply to the vertex coordinate. In one embodiment, the force to apply to a vertex comprises three components, referred to herein as "image force," "torsion force," and "stretch force." The image force tends to cause the vertices to align in relative positions that correspond to "edges" in the three dimensional image produced from the MRA data. Edges in the image correspond to maxima in the gradient magnitude of the image. Roughly speaking, such maxima represent areas where bright pixels are juxtaposed with much dimmer pixels. Such areas are indicative of the luminal surface, e.g. the boundary between the inside of the vessel where blood flows and where the contrast agent responds strongly to magnetic resonance signals, and plaque and tissue where no contrast agent is found. The image force thus tends to cause the vertices of the mesh to assume relative positions that conform to the luminal shape of the vessel.

Production of the isosurface, and application of the image force, may result in a mesh with sharp angles and unevenly spaced vertices. A more useful surface may be obtained by smoothing sharp angles and regularizing the spacing between adjacent vertices. The torsion force may act to improve the smoothness of the mesh. The stretch force may act to regularize the spacing between adjacent vertices in the mesh to create a more uniform appearance.

The manner in which the forces are determined may be better understood with reference to FIGS. 2 and 3. With reference to FIG. 2, a portion 200 of a mesh embodiment comprises four neighboring triangles. Triangles are neighbors when they share at least one vertex. The triangles are EBC, EBA, EAD, and ECD. Each triangle has a normal vector orthogonal to the plane comprising the triangle. The normal of triangle EBC is Ni. The normal of triangle EBA is N . The normal of triangle EAD is N . The normal of triangle ECD is N₂. The vectors Nι-N₄ may be unit normals, that is, vectors having unit magnitude. Vertices are adjacent, e.g. neighbors, when they lie along a same side of a triangle. Thus, vertex E is a neighbor of vertices A, B, C, and D. A set of neighboring triangles comprises a set of neighboring vertices, and may be referred to as a "neighborhood."

The normals of neighboring triangles may be averaged to produce an "average" or "neighborhood" normal for the neighborhood. With reference to FIG. 3, the difference vector N_d between the neighborhood normal N_a and the normal Ni of a triangle 300 is indicative of the sharpness of the angle formed by the triangle 300 on the mesh. The greater the magnitude of N_d, the sharper the angle created by the triangle 300.

A "center" point 302 may be determined for the triangle 300. In one embodiment, the center point 302 has a coordinate with components that are the average of the components of the triangle vertices E, B, and C. In another embodiment, the center point 302 is determined using the well-known center of mass calculation for a triangle (i.e. x = — X_j , where Vi is the set of vertices of

the "' triangle and Γ. is the vector location of the j'^h vertex). In one embodiment, the torsion force to apply to a vertex E is proportional to the magnitude of the "moment arm" vector 308. The torsion force is applied around an axis of rotation 304 through the center point 302.

In one embodiment, determination of the forces to apply to a vertex v follows the diagram of FIG. 4. In general, the image force, torsion force, and stretch force may be determined in other ways. The image force is any force that tends to relocate the vertices of the mesh in a manner such that the positions of the vertices correspond with the vascular surface as identified in the image. The torsion force is any force that tends to smooth sharp edges in the mesh. The stretch force is any force that tends to equalize the distance between neighboring vertices in the mesh. At 412 a gradient magnitude image V7 is derived by convolving the MRA image data I with the gradient of a Gausian function VG :

V/ = VG*7

Selection of the Gausian function is an empirical decision. For example, the non- normalized gradient kernel may take the form:

^2

G ≡ e^~σl

Where r is the distance in centimeters and σ is the space constant of the Gaussian. The space constant may be set, for example, to a dimension (e.g. in centimeters) that is equivalent to 0.5 to 1.5 pixels. Selection of σ is arbitrary, but is selected to be as small as possible to maximize resolution.

→ At 414, a force J on the vertex v is determined by evaluating the gradient of the gradient magnitude image VI at the coordinate of the vertex v: J = [V(V/)](v)

At 416, each unit normal n, of the triangles in the neighborhood of the vertex v are

→ scaled by a value that is proportional to the component of the force J along the

→ → → normal m . The component of the force J along the normal n, is found by taking the dot product:

→ → «,*J

This component is then scaled by an amount c, that is chosen empirically to

produce the best results. The resulting force vector E₍. has the direction of the unit

→ → → normal n, and a magnitude of [c,(«,« J)] :

E_|. =[c,.(«,« J)]«,

→ A force vector F_i is determined for each triangle i in the neighborhood of the vertex

— > → v. At 418 the force vectors F_t are summed to produce the image force vector F, :

_→ neighborhood _→

F, = ∑ F,

To determine the torsion force on a vertex v, an average normal n_a for the

neighborhood of v is determined at 402. The average normal n„ is determined by

summing the normals «, of the triangles in the neighborhood and then dividing by the resulting magnitude to obtain a unit normal vector:

At 404, for a particular triangle i in the vertex neighborhood, a difference

-> vector N_d between the average normal and the triangle normal is determined:

N, «.- «.

At 408 the vector a comprising the axis of rotation for the triangle i is determined by taking the cross product of the triangle normal and the average normal for the neighborhood:

At 406 the moment arm vector / is determined by subtracting the vertex v coordinates from the coordinates of the center o of the triangle I, and then taking the

component of the resulting vector that is perpendicular to the axis of rotation a :

I = p[v -o]

→ Here, p returns the vector component of v - o perpendicular to a . At 410 a magnitude/, of the torsion force to apply to the vertex for the triangle i is

— > determined as the product of the magnitude of the moment arm vector / and the magnitude of the difference vector N_d , scaled by a constant c_t chosen to achieve the best results:

Choice of c, may be made, for example, by matching the shape of the resulting surface to a synthetic phantom, or by matching the surface to the X-ray angiographic image for a vessel of known shape, and is selected to smooth the surface without loss of its shape. Once a value for c_t in a particular situation is determined, it may be used for similar types of scans.

Alternatively, at 420, the torsion force F_t to apply to the vertex v for the triangle i is given the magnitude f_t and a direction determined by the cross product of the axis of rotation and the moment arm, e.g. the direction of the torque is tangential to the arc formed by rotating the moment arm vector around the axis of rotation:

F_t = fΛlx a)

At 422 the torsion force vectors F_t are summed to produce the total torsion force

vector F_τ

_→ neighborho od _→

Fr = ∑ F,

At 424, to determine the total stretch force on a vertex v, a stretch force E. is determined for the vertex v and each neighbor vertex v. . The difference vector between vertex v and the neighbor vertex is scaled by a constant c_s empirically chosen to produce the best results in matching a standard such as a synthetic phantom or a known X-ray angiographic shape:

→ F_s = c_s(v-v_s)

Alternatively, and to minimize shrinkage of the image surface due to the

→ stretch forces, F_s may be defined as the component of the attractive force between vertices that lies in a plane tangent to the surface:

^F, = c, ∑ Proj_P (v - v_s)

J' N

Where N_j ^V is the set of vertices that are adjacent to the j'^h vertex. P. is the plane tangent to the surface at vertex j and proj_p (v - v_) is the orthogonal projection of (v - v.)onto Pj .

At 426, the stretch forces E. between the vertex v and each neighbor are summed to produce the final stretch force for the vertex v:

At 428 the total image force, torsion force, and stretch force are summed to produce the total force on the vertex v:

The total force is determined for each vertex in the mesh, and applied to each vertex to effect relocation of the vertex coordinates. This has the effect of altering the mesh to more model the luminal vascular surface, as indicated by the intensity information in the image, while smoothing the mesh and rendering it more uniform. Vertex relocation may be incremental, in proportion to the magnitude of the total force vector, and in the direction of the total force vector. After all vertices are relocated for a particular iteration, it may be determined that the stopping point has been reached. If so, processing concludes. Otherwise, forces are determined again for another iteration. With reference to FIG. 5, an apparatus embodiment 500 for practicing embodiments of the present invention comprises a processing unit 502 (e.g., a processor, microprocessor, micro-controller, etc.) and machine-readable media 504. Depending on the configuration and application (mobile, desktop, server, etc.), the memory 504 may be volatile (such as RAM), non- volatile (such as ROM, flash memory, etc.) or some combination of the two. By way of example, and not limitation, the machine readable media 504 may comprise volatile and/or nonvolatile media, removable and/or non-removable media, including: RAM, ROM, EEPROM, flash memory or other memory technology, CD-ROM, digital versatile disks (DVD) or other optical storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other medium which can be used to store the desired information to be accessed by the apparatus 500. The machine readable media 504 may be implemented in any method or technology for storage of information such as computer readable instructions, data structures, program modules or other data. Such instructions and data may, when executed by the processor 502, carry out embodiments of methods in accordance with the present invention.

The apparatus 500 may comprise additional storage (removable 506 and/or non-removable 507) such as magnetic or optical disks or tape. The apparatus 500 may further comprise input devices 510 such as a keyboard, pointing device, microphone, etc., and/or output devices 512 such as display, speaker, and printer. The apparatus 500 may also typically include network connections 520 (such as a network adapter) for coupling to other devices, computers, networks, servers, etc. using either wired or wireless signaling media.

The components of the device may be embodied in a distributed computing system. For example, a terminal device may incorporate input and output devices to present only the user interface, whereas processing components of the system are resident elsewhere. Likewise, processing functionality may be distributed across a plurality of processors.

The apparatus may generate and receive machine readable instructions, data structures, program modules or other data in a modulated data signal such as a carrier wave or other transport mechanism. The term "modulated data signal" means a signal that has one or more of its characteristics set or changed in such a manner as to encode information in the signal. This can include both digital, analog, and optical signals. By way of example, and not limitation, communication media includes wired media such as a wired network or direct-wired connection, and wireless media such as acoustic, RF, infrared and other wireless media.

Communications media, including combinations of any of the above, should be understood as within the scope of machine-readable media.

In addition to applications to angiography, the disclosed techniques are applicable to estimation of hemodynamic significance (pressure drops) of stenoses. Such estimations may involve finite element analysis of blood flow, which may depend upon accurate vessel surface models.

Another potential application may be in measuring mechanical stress on the wall of the carotid artery. Such stresses may be related to the development of thromboses (clots). Another potential application that may not involve MRA involves modeling the shape of the colon to detect polyps, where discrimination between a normal colon and a polyp is determined at least in part by the shape of the colon. As used herein, the term "vessel" should be understood to comprise any tubular or approximately tubular surface, not just blood vessels. Examples

Magnetic resonance angiography (MRA) is the magnetic resonance imaging (MRI) of blood vessels in the body. In one type of MRA a contrast agent, such as a gadolinium chelate, is injected into the bloodstream of a patient. The contrast agent causes blood vessels to assert a stronger magnetic resonance signal than the surrounding tissue, and thus to appear brighter in the MRA image. Well known techniques, such as the application of additional magnetic pulses, may be applied to distinguish arteries from veins. Data may be acquired in slices through the region of interest, then combined to form a three dimensional MRA image. Atherosclerotic disease can be detected with magnetic resonance angiography (MRA). As in conventional x-ray angiography, the severity of the disease is determined according to the maximal narrowing of the vessel relative to the normal vessel diameter. Stenosis of the carotid artery, for example, is considered severe if there is a 70% reduction in the vessel diameter in which case surgical intervention is warranted. The degree of stenosis can be measured from the source cross-sectional images, reformatted cross-sectional images (cross-sections at arbitrary orientations) or from maximum intensity projections (MLPs). These visualization methods are all currently available in clinical practice. However, computational measurement of shape, or surface reconstruction, may improve diagnostic interpretation. Firstly, the subjective component in the measurement of vessel shape by visual methods may be reduced using surface reconstruction methods. Vessel surface reconstruction may also be valuable, as a first step, for detection of more subtle vascular disease such as cerebral aneurysms or where vessels are highly overlapped with one another. Finally, vessel surface reconstruction may allow for more comprehensive shape analysis and improved diagnostic accuracy. Mechanical aspects of atherosclerotic disease, for example, can also potentially be assessed by finite-element models constructed from vessel- surface reconstructions and magnetic resonance imaging flow quantification. Such models can potentially be used to estimate wall shear stress and pressure drops at stenoses.

The following examples demonstrate how the disclosed methods may be applied to measure arterial stenoses. These examples are not meant to limit the scope of the invention, both with respect to the methods themselves and their applications. Rather, the examples are presented for illustrative purposes only.

Tables 1 and 2, which appear below, summarize the results for application of the disclosed methods to provide surface reconstructions of arteries and to determine the degree of arterial stenosis. In these Tables, "digital phantom" refers to results for an idealized arterial stenosis blvured to simulate an experimental image, "carotid artery" refers to results for an experimental image of a subject having atherosclerotic disease in the carotid artery, "renal artery 1" refers to results for an experimental image of a subject having noticeable atherosclerotic disease in only the left renal artery, and "renal artery 2" and "renal artery 3" refer, respectively, to the left and right renal arteries of a subject having bilateral atherosclerotic disease in the left and right renal arteries.

1. Synthetic phantom of artery with stenosis A digital phantom of an artery with 70% stenosis was generated in the following manner. The phantom image was generated based on an idealistic axi- symmetric stenosis which is shown in FIG. 6 (d). The slope of the taper of the stenosis relative to the vessel axis was 0.5. The 3D image was generated with isotropic voxel dimensions of 0.12 cm. A high-resolution binary mask of the vessel geometry was then created. The binary mask was then down-sampled 10X in each direction. Each voxel in the resulting image is the sum of the binary values in the original image. This provided an image with the partial volume effects at the edges of the vessel. White image noise was then added to the image. The standard deviation of the white noise is equal to the contrast of the vessel to the background (1000). The image was then blurred by convolution with a Gaussian with a space contant of 1.0 pixel units. A maximum intensity projection of the digital phantom is shown in FIG. 6(a).

The surface of the digital phantom was reconstructed in the following manner. The image was tri-linearly interpolated 2X in each direction to ensure that there was not a cancellation of the image force at the narrowest part of the stenosis from opposing edges. The convolution kernels for computing the gradients for p ^{im ge} _use(j _{a non}__normaij_ze(ι Gaussian function with a space constant of 0.5 pixels. The isosurface was generated for an isointensity value of 2000, and is shown in FIG. 6 (b). Ten thousand iterations were performed. The parameters of the deformable model c_image , c_rotalional , c_stretch used to recreate the digital phantom are shown in Table 1. The units of vertex position were centimeters. The surface reconstruction of the digital phantom is shown in FIG. 6 (c).

The degree of stenosis was manually measured from the 3D display of the deformable isosurface obtained using the disclosed methods. A 64% reduction in the vessel diameter was measured at the narrowest part of the stenosis, representing an error of 6%. In Table 2, the result obtained for the deformable isosurface model is compared with the degree of stenosis determined from non-deformed isosurface image, which showed an error of 14% relative to the true degree of stenosis used to generate the digital phantom.

Table 1. Choice of parameters for each application of the deformable model The mean voxel dimension is the cube root of the voxel volume The scale of the gradient operations refers to σ, the space constant in the non-normalized gradient kernel G

Trial True DegreeTof Isosurface - - Deformable Model Stenosis ^~~ Degree o ^tenosis Degree of Stenosis (error) „ (error!

_LDigitaLPhantom ^≡70% 56% (14%) f ^~ 64% (6%)

_

Carotid-Artery ~^~45% 34% (11%) ^~ 35% (10%)

Renal Artery 1 -14% - 47% (33%) - 8% (6%)

Renal Artery 2 _ 29% 58% (29%) 47%^" (18%)

Table 2. Evaluation of the accuracy of the isosurface and the deformable model according to the degree of stenosis The standard measurements of stenosis (true degree of stenosis) used to calculate errors are (a) the prescribed degree of stenosis for the digital phantom, (b) manual measurement of stenosis from x-ray angiography for the renal arteries and (c) manual measurement of stenosis from the IP for the carotid artery Dependency of the deformable model on the initialization was evaluated as follows Two surface reconstructions of the phantom were obtained using the deformable model with initialization by an isosurface generated a two different isointensity values, 2000 (used for testing discussed previously) and 2500 Silhouettes of the two results (in the orientation shown in FIG. 6) were compared. There was a 3% discrepancy between the two silhouettes where the discrepancy, D, is defined as:

number of discrepant pixels E» = 100 x number of pixels within either silhouette

The discrepancies between the silhouettes were concentrated at the endpoints of the vessel that are not free to deform. At 10,000 iterations, the maximum vertex deformation was 4x10^~6 cm.

2. Renal arteries with atherosclerotic disease

The deformable isosurface was applied to two contrast-enhanced MRA's (Gd-DTPA) of the abdominal aorta and renal arteries with mild stenosis. In the first case (renal artery 1), the deformable model was only applied in the region of the left renal artery since no disease was apparent in the right renal artery (FIG. 7). In the second case, the deformable model was evaluated for both the left and right renal arteries (renal arteries 2 and 3, respectively) since both had evidence of disease (FIG. 8).

The MRA images were acquired on a 1.5T MR system (Philips Medical Systems, Best, Netherlands). The images were acquired with the T1FFΕ sequence. Interpolated voxel dimensions are 0.7812-mm in-plane and 1.9- m slice thickness. Tri-linear interpolation (2x) was applied due to small size of the renal artery. For the first case, the interpolation was applied prior to extraction of the isosurface and computation of the image gradients to maximize the density of vertices for the relatively small renal artery. In the second case, the isosurface was extracted prior to interpolation of the image.

The isosurface was formed by trial-and-error. The primary criteria for selecting the iso-intensity value for the isosurface was that the reconstructed surface had the coπect topology. For these cases, the renal artery is discontinuous at the stenosis if the chosen iso-intensity value is too high. After the isosurface was constructed, the deformable model was only applied to the connected component that constituted the abdominal arteries; background clutter in the isosurface was manually unselected and deleted. This operation required only a single mouse click. The parameters used for the deformable model, c_image , c_rotational , c_stretch are summarized in Table 1 above. 10,000 deformation iterations were performed. The convolution kernels for computing the gradients for E ""^age used a non-normalized Gaussian function with a space constant of 1.0 pixels. The maximum vertex deformations at the 10,000^th iteration were 3xl0^"6 and 7x10^~6 cm for the first and second cases, respectively. The left renal artery, renal artery 1, that is shown in FIG. 7 is mildly stenotic at its origin as confirmed by conventional angiography (FIG 7(e)). The arrow in the MRA image of FIG. 7 (a) indicates the location of the stenosis. As summarized in Table 2 and as shown in FIG. 7, the degree of stenosis of the left renal artery is exaggerated in the isosurface (FIG. 7(b)), whereas the degree of stenosis from the deformable model is more realistic (FIG. 7(d)). The result of surface reconstruction with no torsional smoothing is shown in FIG. 7(c).

The left and right renal arteries (renal artery 2 and renal artery 3, respectively) from a subject with bilateral atherosclerotic disease of the renal arteries are shown in FIG. 8. In FIG. 8(a), the MRA is shown as a sub-volume MIP with arrows indicating the renal artery stenoses. The isosurface initialization is shown in FIG. 8(b), and the result of the deformable isosurface is shown in FIG. 8(c). X-ray angiography, shown in FIG. 8(d), indicated 29% and 48% stenoses of the left and right renal arteries, respectively.

The surface reconstruction by the disclosed deformable model was evaluated according to the correspondence of the apparent degree of stenosis with the degree of stenosis as measured from conventional x-ray angiography. The comparison is summarized in Table 2. The degree of stenoses of the left and right renal arteries as measured from the deformable isosurface was 47% and 35%, respectively. In comparison, the degree of stenoses determined from the isosurface initialization was 58% and 47% for the left and right renal arteries, respectively. 3. Carotid artery with atherosclerotic disease

The disclosed deformable isosurface methods were applied to contrast- enhanced MRA (Gd-DTPA) of a carotid artery with moderate stenosis (FIG. 9). The image was acquired on a 1.5T MR system (GE Medical Systems, Milwaukee,

Wisconsin) using the fast-gradient recalled echo sequence. The image was acquired with a matrix of 256x256 in a 28-cm field of view with 2-ram slice thickness. The image was then zero-fill interpolated 2x in each direction. The isosurface was generated from the image after downsampling the image by a factor of two in each direction since there are not any notable features of the vessel shape that require a mesh with a high vertex density.

The parameters of the deformable model, c_image , c_rotational , c_stretch are summarized in Table 1 above. 10,000 deformation iterations were performed. The convolution kernels for computing the gradients for ρ^ma&e used a non-normalized Gaussian function with a space constant of 1.0 pixels. The maximum vertex deformation at the 10,000^th iteration was 2.2x10^"5 cm.

The surface reconstruction using the deformable model is shown in FIG. 9. Angiographic confirmation of the degree of stenosis was not available. Both the isosurface (FIG. 9(b)) and the deformable isosurface (FIG. 9(c)) are consistent with the MIP (FIG. 9(a)) and with the source images. The deformed isosurface has a smooth appearance. The deformable model noticeably reduces the stair-step appearance produced by the isosurface without changing large-scale features including the degree of stenosis. The degree and relative error of the stenosis determined from the isosurface and the deformed isosurface are shown in Table 2 above.

In view of the many possible embodiments to which the principles of the present invention may be applied, it should be recognized that the detailed embodiments are illustrative only and should not be taken as limiting in scope. Rather, the present invention encompasses all such embodiments as may come within the scope and spirit of the following claims and equivalents thereto.

Claims

I claim:

1. A method comprising: determining image forces for vertices of an isosurface mesh, the image forces determined according to a three dimensional image of a vessel; and applying the image forces to relocate the vertices into relative positions corresponding to a contour of the vessel.

2. The method of claim 1 further comprising: determining torsion forces for the vertices; and applying the torsion forces to relocate the vertices.

3. The method of claim 1 further comprising: determining stretch forces for the vertices; and applying the stretch forces to relocate the vertices.

4. The method of claim, wherein determining image forces for vertices of the mesh further comprises: determining a gradient magnitude image for the image; and for each vertex of the vertices, determining the image force according to the gradient of the gradient magnitude image evaluated at the vertex.

5. The method of claim 2 further comprising: determining, for each vertex of the vertices, the torsion force with a magnitude proportional to the magnitude of a moment vector from the vertex to a center of a triangle comprising the vertex; and determining the torsion force with a direction determined from the cross product of the moment vector and an axis of rotation through the center.

6. A method comprising: selecting an intensity value in an image; producing an isosurface mesh from the image according to the intensity value; and repeatedly relocating vertices of the mesh according to a composite of an image force, torsion force, and stretch force, each determined for each of the vertices, until a stopping point is reached.

7. The method of claim 6 further comprising: reaching the stopping point when relocating the vertices fails to produce substantial changes in the positions of the vertices.

8. The method of claim 6, wherein determining the image force for each of vertices further comprises: determining a gradient magnitude image for the image; and for each vertex of the vertices, determining the image force according to the gradient of the gradient magnitude image evaluated at the vertex.

9. The method of claim 6 further comprising: determining, for each vertex of the vertices, the torsion force with a magnitude proportional to the magnitude of a moment vector from the vertex to a center of a triangle comprising the vertex; and determining the torsion force with a direction determined from the cross product of the moment vector and an axis of rotation through the center.

10. The method of claim 6, wherein the stretch force for each vertex of the vertices is proportional to distances between the vertex and its neighboring vertices.

11. The method of claim 6, wherein the stretch force for each vertex is proportional to a component of the displacement between vertices that is tangent to a surface containing the vertices.

12. A method for reconstructing a shape of a vessel in a subject from a magnetic resonance image, comprising: selecting an intensity value in the image; producing an isosurface mesh from the image according to the intensity value; and repeatedly relocating vertices of the mesh according to a composite of an image force, torsion force, and stretch force, each determined for each of the vertices, until a stopping point is reached to provide a reconstructed shape of the vessel.

13. The method of claim 12 further comprising detecting a stenosis in the vessel from the reconstructed shape of the vessel.

14. The method of claim 12 further comprising: reaching the stopping point when relocating the vertices fails to produce substantial changes in the positions of the vertices.

15. The method of claim 12, wherein determining the image force for each of vertices further comprises: determining a gradient magnitude image for the image; and for each vertex of the vertices, determining the image force according to the gradient of the gradient magnitude image evaluated at the vertex.

16. The method of claim 12 further comprising: determining, for each vertex of the vertices, the torsion force with a magnitude proportional to the magnitude of a moment vector from the vertex to a center of a triangle comprising the vertex; and determining the torsion force with a direction determined from the cross product of the moment vector and an axis of rotation through the center.

17. The method of claim 12, wherein the stretch force for each vertex of the vertices is proportional to distances between the vertex and its neighboring vertices.

19. The method of claim 12, wherein the stretch force for each vertex is proportional to a component of the displacement between vertices that is tangent to a surface containing the vertices.

20. The method of claim 12, wherein the vessel is a colon and further comprising detecting a polyp in the colon from the reconstructed image.