WO2000052629A1

WO2000052629A1 - 3-d shape measurements using statistical curvature analysis

Info

Publication number: WO2000052629A1
Application number: PCT/US2000/005596
Authority: WO
Inventors: John E. Stewart
Original assignee: Virginia Commonwealth University
Priority date: 1999-03-03
Filing date: 2000-03-03
Publication date: 2000-09-08
Also published as: CN1399763A; CA2364176A1; AU3720700A; JP2003502723A; EP1236156A1

Abstract

A three dimensional curvature algorithm (fig. 1) using linear regression for modeling biological matter. A three dimensional representation of a structure is obtained by scanning the matter. Selected regions of the scanned structure are assigned numerical values. A number of curvature measurements are calculated based on vertices. Linear regression analysis is used to obtain a coefficient of regression for all curvatures. Variance inflation factors are calculated for the curvature measurements. Multiple regressions are performed to obtain a best fit model.

Description

3-D SHAPE MEASUREMENTS USING STATISTICAL CURVATURE ANALYSIS

DESCRIPTION

BACKGROUND OF THE INVENTION

Field of the Invention

The present invention generally relates to a computer implemented method for applying differential computational geometry and statistics to detect shapes on three dimensional (3-D) computer models representative of a sequential series of two dimensional images of a structure. Twelve new curvature measures are also described. These measurements are applied to the surface of a 3-D computer model made up of triangles. Collections of 3-D images of known structures having aberrant surface curvature characteristics associated with dysfunction or defects are then used to train the software to recognize the particular curvature properties of an aberrant structure. Multiple linear regression is applied to the curvature measures of the 3-D test models and the results are optimized to remove multicolinearity and maintain a high coefficient of regression.

Background Description

An integral step in many medical diagnostic protocols is the detection of structural abnormalities associated with anatomical regions of interest.

Accordingly, much of present day medical diagnosis involves the structural analysis of medical images either directly, through interpretation of film images, or indirectly, through the generation and interpretation of 3-D computer models of an anatomical area of interest. Analysis of medical images can be difficult and time consuming, with the subjective view of the observer often clouding the interpretation of such images. A computerized method for providing objective measurements of the surface area and shape of an internal anatomical structure would therefore provide improved detection of abnormalities associated with structural defects.

By way of example, cerebral aneurysms, also known as "berry aneurysms," are typically found as sacular outpoutchings of major intracranial arteries. They are known to enlarge progressively in many individuals and are thought to arise from arterial blood flow striking a weak area of the internal elastic lamina in the wall of an artery. Based on autopsy studies, these aneurysms occur in the general population at a rate anywhere between 1-8 %. For a given individual with an unruptured aneurysm, the annual risk of rupture is 1-2% while the lifetime risk of such a rupture is approximately 50%. At present time there do exist methods for elective treatment of such aneurysms in order to prevent rupture. Accordingly, effective diagnosis can potentially greatly reduce the morbidity and mortality associated with the occurrence of aneurysms.

Traditional diagnostic protocols for detecting an aneurysms include cerebral angiography and magnetic resonance angiography (MRA). Cerebral angiography is an invasive procedure in which a catheter is advanced through the intracranial vessels and a contrast material is injected to enable visualization of the vessel walls. This procedure is costly, time consuming, and involves serious potential complications that can lead to aneurysm rupture, stroke, and death, although the risk of such complications is small (less than 1%). MRA, on the other hand, is less invasive, less time consuming, and presents less risk than cerebral angiography.

Differential geometry has been used extensively as a means of analyzing a variety of geometric shapes and computational geometry is a means of rapidly determining surface characteristics using the mathematics of differential geometry. The techniques of computational geometry have been applied to the prediction of protein-protein interactions (Duncan et al, "Shape Analysis of Molecular Surfaces," Biopolymers, 33: 231-238, 1993), analysis of dose-effect surfaces of combined agents (Lam, "The Combined Actions of Agents Using Differential Geometry of Dose-effect Surfaces," Bulletin of Mathematical Biology, 54:813-826, 1992), and the study of biological surface growth.

Computational geometry has been applied in the medical field in the context of curvature analysis to diagnose lung cancer (Kawata et al., "Computer Aided Diagnosis System for Lung Cancer Based on Helical CT Images," ICIAP '97 Proceedings, 2: 420-427, 1997), cardiac aneurysms

(Fantini et al., "Quantitative Evaluation of Left Ventricular Shape in Anterior Aneurysms," Catheterization and Cardiovascular Diagnosis, 28: 295-300, 1993), and assessment of the treatment of brain tumors (Dai et al., "Intracranial Deformation Caused by Brain Tumors: Assessment of 3-D Surface by Magnetic Resonance Imaging," IEEE Transactions on Medical

Imaging, 12, 4:693-702, 1993). All of these procedures have employed smoothing in order to avoid the large variations in curvature that are often seen when measuring the curvature of anatomic structures. The smoothing process, however, often distorts the measurement of curvatures. Recent developments have applied surface curvature analysis to extract feature lines on 3-D surfaces (Monga et al, "From Partial Derivatives of 3-D Density Images to Ridge Lines," Visualization in Biomedical Computing, VBC'92, Proceedings, pp. 118-129, 1992) but not as applied to 3-D surfaces of the cerebral vasculature. MRA studies of the brain commonly consist of 72-90 two dimensional

(2-D) cross sectional images. Arteries in these images are highlighted using a 3-D time-of-flight algorithm or its equivalent that is included as a matter of course with typical magnetic resonance scanner software. The images are displayed on sheets of 2-D radiological film for review on light boxes. Each image is read by a radiologist or surgeon who "creates" a three-dimensional interpretation of the vasculature through his or her knowledge of the anatomy of the blood vessels. This image may not be reduced to a three dimensional interpretation that can be reviewed by third parties. Due to the branching nature of the vasculature, this interpretation can be very difficult and time consuming. Branching vessels often appear as an aneurysm on a single image but analysis of surrounding images reveals normal vasculature. Similarly, an aneurysm might appear as a branching vessel unless it is carefully followed through a series of images to its termination.

The use of computers to assist in detecting cerebral aneurysms prior to this invention has been limited. Attempts have been made to analyze image data to detect cross-sectional images representative of aneurysms but these attempts have been futile, due to the difficulty in distinguishing an aneurysm from a normal branching point for an artery (Fessler et al., "Object-Based 3-D Reconstruction of Arterial Trees from Magnetic Resonance Angiograms," IEEE Transactions on Medical Imaging, 10:1, 1991). Also, skeletal representations of the vasculature have been used to determine vessel branch points and topological features in vessels (Puig et al., "An Interactive Cerebral Blood Vessel Exploration System" Visualization '97, Proceedings, pp.443- 446, 1996). While these techniques are useful in some fields, such as virtual colonoscopy, the small vessel radii and large directional changes encountered in the topography of blood vessels are problematic to the application of skeletal representations for the identification of aneurysms. Smoothing and filtering of MRA images have been employed in an attempt to overcome this problem but the result is often a distorted view of the vasculature. Morphometric analysis has been proposed as a means of automatically detecting aneurysms but this procedure can only analyze small regions of the vasculature at a time (Matsutani, et al, Quantitative Vascular Shape Analysis for 3-D MR-Angiography Using Mathematical Morphology," Computer Vision, Virtual Reality and Robotics in Medicine, CVRMed '95, Proceedings, pp.449-454, 1995).

SUMMARY OF THE INVENTION

Based on the foregoing, a novel method for generating a more objective 3-D representation of anomalous regions on the surface of structures is disclosed. The present invention provides a method for objectively determining a 3-D analysis of a series of 2-D images. One particular embodiment of the invention is based on the simple observation that cerebral vessels are roughly cylindrical while aneurysms contained within such vessels are roughly spherical in shape. This particular embodiment exploits the spherical nature of aneurysms and computational differential geometry to identify and highlight aneurysms in unread 2-D MRA images.

There are two classical mathematical entities used to analyze smooth surfaces which are known as the first and second fundamental forms. These measures are useful because they are intrinsic to the surface and therefore invariant to transform (rotation, translation, scaling). The mean and Gaussian curvatures are based upon the fundamental forms.

In classical differential geometry, the principal curvatures k, and k₂, are averaged to produce the mean curvature (H) and multiplied to produce the Gaussian curvature (G). The mean and Gaussian curvatures can be used to classify the surface into particular types. In general, if H>0, a surface is concave and if H<0, the surface is convex. Beyond this, surface classification can be broken down into elliptic, hyperbolic, parabolic or planar according to the following rules. If G>0, the surface is elliptic, if GO, it is hyperbolic, if G=0 but H ≠ 0, the surface is parabolic and if G=0 and H=0, the surface is planar. It should be noted that cylindrical surfaces are parabolic and spherical surfaces are elliptic. Most algorithms that classify real surfaces according to curvature utilize an approach known as thresholding, wherein the principal curvatures are represented as a continuum and a particular threshold is used to separate the surface into different curvature types.

One particular technique has been utilized previously to approximate the principal, mean and Gaussian curvatures at the vertices of a triangulated surface. This technique is based on the fact that a surface can be locally represented as a graph of a bivariate function and the curvature of this bivariate function can be directly measured (Hamann et al., "Curvature Approximation for Triangulated Surfaces," Computer Suppl. 8: Geometric Modeling, Springer-Verlag Wein, New York, pp. 139-153, 1993). If the function provides a good approximation of the surface shape, then the curvature measured on the bivariate function will also be a good approximation of the underlying surface curvature. Since the interest in this process is the curvature at each vertex, a quadratic polynomial is fit to each platelet (region surrounding a vertex).

The method of measuring curvature described previously (Hamman et al.) allows for a platelet to consist of only a single row of vertices around a central vertex (Figure 2). This single row of vertices is then used to determine the bivariate quadratic polynomial at each surface vertex. On "ideal surfaces" that are very smooth, this approximation works well. For surfaces that are not ideal, ones that have bumps and triangles of poor aspect ratio, this single row of surrounding vertices does not provide adequate data to determine a bivariate quadratic polynomial that accurate approximates the local surface shape. Therefore, the technique described above has been modified to allow multiple rows of surface vertices to be added to the platelet and used in the determination of the bivariate quadratic polynomial. The result is a much more accurate approximation of local surface shape and subsequently a more accurate approximation of surface curvature around each vertex. The determination of the number of additional rows of surface vertices to include is optimized by finding the number of rows of surface vertices that produces the largest coefficient of regression (described later). The classical curvature approximations provide a great deal of information about the surface shape and rate of change in surface shape. The classic measures are not, however, optimized to recognize and distinguish between cylinders and spheres. Surface size will also confound the application of classical curvature methods to cylinders and spheres. For example, a large spherical object may have the same Gaussian curvature as a small cylindrical object. Similarly, large and small objects of the same shape will have dramatically different curvatures, making it difficult to locate aneurysms, for example, with these methods independent of their radii.

Twelve new curvature measures are disclosed that are useful to distinguish between spherical and cylindrical images presented in a 2-D format.

Furthermore, these novel curvature measures are normalized to allow spheres to be recognized independent of their radii. Optimization of these measures through the application of multiple linear regression provides an optimum subset of curvatures that can be used to extrapolate structural anomalies from triangulated surfaces taken from 2-D image data sets.

It is a further object of the invention to provide a method of evaluating determining three dimensional structures comprising the steps of: a) obtaining a computerized three dimensional representation of a structure or structures; b) identifying a first set of regions on the three dimensional structure or structures and assigning a numeric value to said structure or structures; c) identifying a second set of regions and assigning a numerical value to said regions; d) determining values for a plurality of curvature measures for each vertex on a surface of the structure or structures; e) performing multiple linear regression analysis on said values determined in said determining step to obtain a coefficient of regression for all curvatures for all vertices; f) determining the variance inflation factor for each of said curvature measures; g) if all variance inflation factors are less than 10, go to step 1; h) if any variance inflation factor is greater than 10, sequentially reduce the subset of curvature measures used in multiple linear regression by 1 ; i) performing multiple linear regression on all combinations of curvature measures possible for each subset; j) selecting the subset of curvature measures that yields the largest coefficient of regression; k) performing multiple linear regression analysis on said values determined in d) to obtain a coefficient of regression for said curvature subset; and 1) inserting the partial coefficients of linear regression into the linear equation generated by multiple linear regression for said curvature subset.

BRIEF DESCRIPTION OF THE DRAWINGS

The foregoing and other objects, aspects and advantages will be better understood from the following detailed description of a preferred embodiment of the invention with reference to the drawings, in which:

Figure 1 is a flow chart representing a method in accordance with the present invention of applying novel and classical curvature measures and multiple linear regression to optimize structural differences in a series of sequential vertices. Figure 2 is a perspective view of the first principal curvature vector.

Figure 3 is a perspective view indicating the normal vectors that comprise a normal triangle on the surface of a unit sphere.

Figure 4 is a perspective view indicating surface normal vectors lying in a plane perpendicular to a cylinder. Figure 5 is a perspective view indicating the normal triangle characteristics for a sphere.

Figure 6 is a perspective view showing the radius of an inscribed and a circumscribed circle in relation to the normal triangle.

Figure 7 is a graphical representation of the unit normal vectors of a surface triangle and normal triangle.

Figure 8 is a graphical representation of the correlation between the correlation coefficient and platelet radius.

Figure 9 is a perspective view of a platelet constructed from vertices and triangles.

DETAILED DESCRIPTION OF A PREFERRED EMBODIMENT OF THE INVENTION

The present invention generally relates to a method and system as schematically represented in Figure 1 , of applying differential computational geometry to the analysis of surface curvature. As shown in Figure 1, the method accepts as input a 3-D computer model of a structure that can be generated either automatically or semi-automatically from a collection of cross sectional images.

Medical applications of the invention disclosed herein can be adapted to the study of much of the human body, as well as the bodies of other animals. Much of the human body is composed of organs that are roughly spherical or roughly cylindrical in shape. In many instances, derivation from one of these shapes is strongly correlated with disease. For example, the present techniques can be applied to the measurement of arterial stenosis, the detection of arteriovenous malformations, colonic polyps, and the detection of lung and liver cancers.

The present invention can also be applied to the analysis of spherical and cylindrical machine parts that cannot otherwise be readily taken apart for inspection. In this particular embodiment, non-destructive testing of machinery via detecting anomalies in structure of internal parts of a machine. The invention is also applicable to the analysis of 3-D imagery generated from weather patterns to detect funnel clouds or thunderstorms, for example. In another embodiment of the invention, the method can be applied to the technique of molecular modeling to detect abnormal protein-ligand interactions and in confocal microscopy detecting cancerous cells in a collection of images taken from a microscopic field or fields.,

Sequential cross sectional images of a selected anatomical structure are acquired with the use, for example, of a scanner such as a magnetic resonance imaging scanner (MRI). Any type of digital image scanner can be used in place of MRI images such as a helical computer tomography scanner (CT), ultrasound, or PET images. The 2-D images are arranged in a computer memory to create a 3-D data volume set. The image data to be analyzed can be generated and or stored in any of a variety of image formats. For example, the present invention is ideally suited for use with Picture Archiving and Communication system (PACS) format. The image data can be stored in the digital imaging and communications in medicine standard (DICOM), or as raw binary slices, or in a variety of volume formats. The image data can be stored in the computer memory in an internal data format which allows the image files to be saved as a single data volume instead of individual image files. The internal data format can also permit standard compression and uncompression techniques to be incorporated, reducing computer disk storage needs

In one embodiment of the invention, a Web-based PACS software package acts as an automated filing system accepts and stores digital images created by traditional 2-D means. A second software component (IsoView) reads DICOM images stored in the PACS format and creates manifold (closely and singly connected )3-D triangulated surfaces from these images using a variant of the marching cubes algorithm. The surface images can be interactively displayed on a local workstation, stored as Virtual Reality Modeling Language (VRML) for later review, or stored as video recording or photographs for future viewing. The 3-D images created in the IsoView format serve as the input for the curvature measurement method disclosed herein.

The calculation of the bivariate quadratic polynomial at each surface vertex is computed exactly as described by Hamman et al. The principal curvature measures, kl and k2, are then computed on these bivariate quadratic polynomials as described in the same. Hamman's technique determines the principal curvature magnitudes, kl and k2, but does not determine the principal curvature vectors, vl and v2. Therefore, equations and a mathematical method are described below to determine the principal curvature vectors vl and v2. The first step is to compute the eigenvectors (which were not described or computed by Hamman, et al.) of the same matrix for which the eigenvalues were determined,

Equation 1

'2,0

- A = '1

'U ^■O,2

The eigenvectors are calculated by solving,

quation 2

The solution to this system of equations is,

Equation 3

The 3-D principal curvature vectors vl and v2 are computed by normalizing the parametric principal curvature vectors and multiplying them times the 3D basis vectors determined in the method as described by Hamman,

V_j = e_nιb₁ +e_Vι b_{2 5} j = l, 2 Equation 4

The result are the two principal curvature vectors, vl and v2, that are used in newly disclosed curvature measures 3, 4, and 5 discussed below.

As noted above, the classical curvature approximations provide a great deal of information about the surface shape and rate of change of surface shape, they are not optimized to recognize and distinguish cylinders from spheres. In addition, curvature magnitude also masks surface curvature analysis with the classical approximations. The present invention describes twelve novel scalar measures that are applicable to measuring curvatures. Two of these scalars combine the Gaussian and mean curvatures into a single value, three measure the change in principal curvature vector direction and seven make use of the surface normal vectors to predict surface shape.

The calculation of these novel curvature measures are described as follows, where "||" represents the Euclidean norm:

1. MIN(||kl ||,||k2 ||)/MAX(]|Kl,||k2 ||): A ratio of the two principal curvatures is calculated. Since both kl and k2 can assume positive (concave) or negative (convex) values, the unsigned magnitude of kl and k2 is determined. The smaller of the two numbers becomes the numerator and the larger becomes the denominator. This ration should be 0 for a cylinder, 1 for a sphere, and undefined for a plane.

2. ||kl-k2|| : The difference between the two principal curvatures is calculated. This value should be 0 for spheres and planes and 0 otherwise. Referring now to Figure 2,

3. AVGdlV_jβVi_j:!), j=l,n: The average of the dot product between the first principal curvature vector 200 at the center platelet vertex 201 and the principal curvature vectors 202 of those vertices immediately connected to this vertex is calculated. For convex surfaces, 200 will be parallel to the axia of a cylinder and will be randomly oriented on the surface of a sphere. Therefore this curvature measure should be small for cylinders and large for spheres. It also will be large in regions where the surface shape is changing dramatically such as at the branch point of vessels.

4. AVG(||v₂*v_2j:||), j=l,n: The average of the dot product between the second principal curvature vector at the center platelet vertex 201 and those vertices immediately connected to this vertex is calculated. For convex surfaces, v₂ will be perpendicular to the axia of a cylinder and will be randomly oriented on the surface of a sphere. Therefore this curvature measure should be large for spheres, small for large cylinders, and large for small cylinders. It should also be large at vessel branch points.

5. AVGfllV_jΦVu.l), j=l,n) - AVG(||v₂»v_2j:||), j=l,n: The difference between the two curvatures in 3 and 4 above is taken. This measure should be small for small cylinders and large for large cylinders and spheres.

Seven more curvatures disclosed herein use the surface normal vectors to predict surface shape. Referring now to Figure 3, all of these measures rely on one simple observation, i.e., for surface 300 as shown in Figure 3, if the normal vectors n₀ n, n₂at the three corners of a triangle 301 which lies on surface 300 are drawn such that they originate from a single point, these normal vectors will form a triangle . For unit normal vectors ri_j, the vertices of this new triangle will lie on the surface of unit sphere 302. The circumference, area and aspect ratio of this triangle 303 ( the "normal triangle") can provide considerable insight into the shape of the surface from which the triangle originated. Figure 5 illustrates the characteristic of the normal triangle for a cylindrical surface. A cylinder is defined by a circle extruded along a vector running perpendicular to the plane of the circle. Referring now to Figure 4, surface normals 401, 402, and 403 for a cylinder lie in planes perpendicular to the axis of the cylinder. If the unit normal vectors for three points on the surface of a cylinder are translated towards one another parallel to the axis of this cylinder, the three planes will eventually coincide and all three unit normal vectors will lie in this coincident plane. The "tips" of these unit normal vectors will form a circle 404 with a radius 1 unit larger than the radius of the cylinder 405. If each unit normal vector is then translated in 2-D with the coincident plane, they can be positioned such that they originated from the same point but still existed entirely within the coincident plane. The tips of these unit normal vectors will therefore form triangle 406 that lies in the coincident plane.

This normal triangle will typically have a large aspect ratio that will increase to infinity as the dot product of the normal vectors approaches 1.0. In other words, the more triangles on the surface of the cylinder, the larger the aspect ratio of the cylinder's normal triangles. Because the radius of the cylinder does not influence the aspect ratio of the normal triangle, this measure will be useful in detecting cylinders of any radius.

Referring now to Figure 5, which illustrates the normal triangle characteristics for a sphere. All surface normals for a sphere are directed in a path from the center of the sphere through a surface vertex. If a plane is formed such that any two surface vertices and the sphere center are members of this plane, then the two unit normal vectors at the surface of the sphere will lie within this plane. The center of the sphere, the two surface points, and the tips of the two unit normal vectors will form two triangles 501 and 502, as shown in Figure 5. Equation 5, x_n=(x_t(r+l))/r defines, using the similar triangles theorem, the length of x_m - one side of the normal triangle. If the same steps are carried out for the two remaining sides of the normal triangle, it becomes clear that each side is increased in proportion to (r+l)/r from the length of the surface triangle. If all three sides of a triangle increase by the same proportion, the triangle formed (normal triangle) will have the same aspect ratio as the original surface triangle. Therefore the aspect ratio of the normal triangle depends on the aspect ratio of the surface triangle and is independent of the radius of the sphere. A better way to make use of the normal triangle curvature is to construct an equilateral triangle around the surface vertex of interest and find the normals at the vertices of this new triangle. Because the bivariate quadratic polynomial produces a good approximation of surface shape an equilateral triangle is created on this polynomial surface. Therefore the surface triangle area is uniform and does not influence normal triangle area. The normals at the vertices of this equilateral triangle are calculated using the quadratic polynomial coefficients shown below in Equation 6:

r_u

The size, aspect ratio, and tilt of the normal triangle are used to predict the shape of the underlying triangles. For a spherical surface, the normal triangle aspect ratio will be the same as the surface triangle aspect ratio and will therefore also be an equilateral triangle, as shown in Figure 6. For a cylindrical surface, the normals will fan out in a straight line forming a triangle orthogonal to the axis of the cylinder. If the surface under examination is planar, all three normals will coincide and their tips will form a point. The size of the normal triangle also reflects the magnitude of curvature on the surface. Large normal triangles are seen in highly curved regions while small normal triangles are seen in relatively flat regions. Curvature measures that use these normal triangles are referred to as "normal triangle curvatures."

The seven additional novel scalar measures are:

6. Radius of an Inscribed Circle (inside the normal triangle). (INSC) Referring now to Figure 6, there is only one way to fit a circle 602 inside the normal triangle 600 such that the three triangle sides are tangent to the circle and the circle center is at the intersection of the three angle bisectors. The three sides of the triangle have lengths L₀, L„ and L₂The radius of this circle 601 is calculated and should be large for a sphere, small for a cylinder, and 0 for a plane:

S = -(L₀ + L_x + L₂ ) Equation 7

(S-L₀)(S-L.)(S-L₂) Equation s t ^'

S

7. Radius of a Circumscribed Circle(CIRC):Referring again to Figure 6, there is also only one way to fit a circle 603 through each of the normal triangle vertices such that the center is at the intersection of the perpendicular bisectors of the three sides of the normal triangle 600. The radius 604 of this circle 603 is calculated and should be very large for cylinders, 0 for planes, and small for spheres:

R_ci_rc ⁼ - — ^" — ^~ Equation 9

4 R_tasc S

8. Radius of an Inscribed Circle/Radius of a Circumscribed Circle (RINCIR): This ratio measures the aspect ratio of the triangle. It should equal 0.5 for a perfect sphere and nearly 0 for a cylinder. This ratio is undefined for a plane. This ratio is relatively insensitive to curvature magnitude since this magnitude will appear in both th enumerator and denominator. It should therefore be large for spheres and small for cylinders independent of their respective radii.

9. Area (NAREA): The area of the normal triangle is calculated by taking the cross product of two sides of this triangle and dividing by 2.0. It should be large for a sphere, small for a cylinder and 0 for a plane.

10. Perimeter (NPER): The perimeter of a normal triangle should be large for a sphere, of intermediate length for a cylinder and 0 for a plane.

11. Area/Perimeter² (RNAP): Like measure 8, this measure is used to find the aspect ratio of the triangle. It is large for a sphere, small for a cylinder and undefined for a plane. Because the perimeter does not vary in magnitude as much as the circumscribed circle, this measure is more stable on cylindrical surfaces than measure 8.

12. TILT: TILT is measured by first finding the unit normal vector of the surface triangle. This can be done by calculating the cross product of two of the three sides of the surface triangle and then normalizing (converting the vector to unit length) the normal vector of the surface triangle. The unit normal vector of the surface triangle is illustrated in

Figure 7 as n_s. Similarly, the unit normal vector, n_nt of the normal triangle is also determined. TILT is calculated by finding the dot product of these two vectors, n_s and n_nt This scalar value will be 1 for convex spherical surfaces, 0 for cylindrical surfaces and -1 for concave spherical surfaces.

The above equations produce a group of statistical measures to determine curvature. A useful technique for combining multiple variables into a single equation is multiple linear regression. Multiple linear regression minimizes the sum of the squares of the residual between the regression equation and data to produce the linear equation that best fits the data.

Furthermore, multiple linear regression uses the sum points over the entire surface. Significant variations in surface curvature ("Noise") seen in some curvature measures can be averaged out in this process. All that is required is for an underlying pattern to exist for the surface as a whole. Another attractive feature of this technique is that the predictive ability of the equation improves as more 3-D surfaces are analyzed. In other words, the aneurysm predictions will become more accurate as more aneurysms are entered into the equation.

Multiple linear regression is used to fit two or more variables to a set of data. The objective of the technique is to determine the coefficients to y=j3₀ + j3,X, + j3₂X₂ + ...β_mX_m +e, Equation 10 that produce the minimum sum of the squares of the residuals around the regression line as,

Equation 11

In these equations, y is the dependent variable, x_m are the independent variables and β are the partial coefficients of linear regression. Differentiating the above equation with respect to each of the regression coefficients results in the ma trix

Equation 12

The x_mι are the 16 curvature approximations consisting of the four classical curvature measures and the twelve novel measures reported herein, while y, are values assigned to the vertices of the 3-D computer models that "teach"multiple linear regression what aneurysms look like.

EXAMPLE 1 In order to test the usefulness of the disclosed techniques in predicting aneurysms, MRA studies of 11 patients with normal arterial vasculature and 11 patients with diagnosed aneurysms were reconstructed. The 3-D models of all 22 patients were created at a grayscale threshold of 325. The total number of vertices for all 3-D computer models is 290,802. The aneurysm patient group was contained nine females and two males while the normal group consisted of eight females and two males. The average age of the aneurysms group was

53.2±17.7 years with a range from 22 to 78 years. The average age of the normal group was 44.6±15.3 years with a range of 24-73 years.

After reconstruction, the surface vertices belonging to the aneurysm group were partitioned into three groups- aneurysm, transition region, or normal. This was done interactively with IsoView. The entire dome of each aneurysm was partitioned and given a value of 1.0. The transition regions for each aneurysm were marked separately and given a value of 0.5. All other points on the model were assigned a value of 0.0. All vertices that make up the 3-D computer models of the normal patients were assigned a value of 0.0. These values constitute the y, dependent variables of the two previous equations. The transition regions lying between the aneurysm and normal vasculature is included because it has different surface curvature characteristics than the aneurysm dome and is often seen accompanying an aneurysm. The transition region is typically concave and is commonly referred to as the neck of the aneurysm.

Having defined all variables, multiple linear regression of all data is then performed. The first step in this process is to compute all 16 scalar curvatures at each vertex of all 22 3-D computer models. The curvature approximations of equations 1-12 are measured for each model and stored in separate files for each patient and curvature type. A platelet radius determines the number of points to use in the least squares fit. The larger this radius, the more points are used. This modification solves the problem of directional bias that can occur when the image pixel size is much smaller than the distance between consecutive MRA images. Directional bias would cause the preferential selection of points in the direction normal to the image planes since triangles would tend to be "stretched" in this direction. No matter what radius is chosen, the first concentric row of points around the central vertex are always used.

Because the regression of 16 variables over 290,802 vertices may be too large a problem to solve utilizing standard statistical software, software was developed and optimized specifically to solve this large multiple linear regression problem. Such software was validated against previously solved problems found in Freund et al, Statistical Methods. Academic Press Ltd., San Diego, Ca, 1997 and Chapra et al, Numerical Methods for Engineers with Personal Computer Applications. McGraw-Hill, Inc, New York, New York, 1985.

A useful measure of the goodness fit of the regression line to the data is the correlation coefficient as in,

Equation 13

where Sr is given by Equation 11 ,

^s

-øo -fcrfe, --ftA,⁾²

and S_t is the sum of the squares around the mean for the dependent variable yj. The question then becomes what platelet radius to use in measuring the curvature on the 3-D computer models. The simplest way to answer this question is to try different radii to determine the radius that results in the largest correlation coefficient.

The 16 scalar curvature measures for all 290,802 vertices were computed for platelet radii of 0.0 mm (i.e., one concentric row of points) to 8.0 mm in increments of 1.0 mm. Multiple linear regression was then performed on each group of curvatures and the regression coefficient was computed. Referring now to Figure 8, the relationship between correlation coefficient and platelet radius for the multiple linear regression of all 22 3-D computer models using all 16 scalar measures of curvature is shown. Clearly a platelet radius of 3.0 mm produces the best correlation coefficient: 0.17011. Although this correlation coefficient does not indicate a perfect linear relationship, a test of the model was performed. The test for the model is simply a test of the hypothesis that the entire set of coefficients associated with the 16 independent variables of equation (5.1) is 0. The alternative to this hypothesis is that one or more of the coefficients are not 0. The F statisticis defined by,

is used to evaluate this hypothesis where r2 is the coefficient of determination (equal to the correlation coefficient to the power of 2), m is the number of curvature measures and n is the number of vertices analyzed. The null hypothesis is,

Ho: (βi, β₂, .... β_m)= 0 Equation 15

The F statistic for our regression model using all vertices and all curvatures at a platelet radius of 3.0 mm is, F(16,290785)=541.55.

The p=0.005 critical value for F(16,∞) is 2.14. Therefore, the null hypothesis that the regression model does not exist can be rejected and the partial coefficients β for the independent variables x are not all 0. The next step in confirming the model is a test of the individual partial regression coefficients. The test statistic for this hypothesis is also the F statistic but the equation to test a single partial coefficient 3. is,

F _EY(il,n -m - _Il_\) A = —H - — J —J I Equation 16

MSE

where c_}J is an element along the diagonal of the inverted multiple regression m x m matrix of equation (5.3) and MSE is the mean standard error defined by, Equation 17

1VISK n — m

The null hypothesis is,

H₀:j3.O.

The p=0.005 critical value for F(l,∞) is 7.88. Table 1 lists the F statistic, the null hypothesis test and the partial ( .) and standardized (β *) correlation coefficient for the multiple regression of all models and all curvatures for a platelet radius of 3.0 mm. The standardized partial correlation coefficients are defined by s β * -=β -----A . Equation 18

where s_x, is the sample standard deviation of xi and sy is the standard deviation of y. These coefficients provide a means of comparing coefficients across variables by standardizing the means to 0 and variances to 1.

TABLE I. Hypothesis testing and partial regression coefficients for each curvature type.

The above table discloses that the curvatures G, CIRC, RINCIR, NPER, RNAP, and TILT should contribute the most to the prediction of aneurysms. These six equations passed the null hypothesis. The preceding analyses demonstrate that a subset of curvature measures are useful in predicting aneurysms but does not provide any information as to how these variables are related to each other. Multiple linear regression attempts to determine the change in the dependent variable for a given change in an independent variable while holding all other independent variables constant. If two independent variables are colinear, it will be impossible to hold one constant while changing the other. The result of this violation is that one or more of the partial regression coefficients may attempt to describe a phenomena that is not exhibited by the data. This phenomenon is known as multicolinearity.

A useful measure of multicolinearity is the variance inflation factor (VIF), which indicates how much larger the variance of each independent variable is than it would be if it were not correlated with other independent variables. VIF is calculated by first performing a multiple linear correlation between each independent variable and all other independent variables. One of the independent variables x. is removed from the left hand side of equation 5 and is used in place of y. The coefficient of determination (r j ) ² is then determined for each independent variable and used to compute VIF:

VIF = — ¹ 2 ^rj Equation 19

A VIF of 1 indicates that independent variable x. is not multicolinear with other independent variables. Anything larger than 1 indicates some degree of multicolinearity. Although there are no firm rules as to the magnitude that VIF can attain before multicolinearity is demonstrated, many statisticians have adopted a value of 10 as the VIF cutoff. Thus, independent variables greater than 10 are considered multicolinear while those less than 10 are considered independent.

Table II lists the VIFs for the five classical curvature measures and the twelve novel measures of curvature described herein. Columns labeled subset size 16-5 represent the test of the correlation coefficient for a combination of

16-5 combinations of curvature subsets. The curvatures were calculated on a bivariate quadratic patch created with a platelet radius of 3 .0 mm. Aside from CIRC and TILT, all curvatures are strongly multicolinear with each other when all 16 curvatures are analyzed. Having demonstrated a regression model, it must be determined which independent variables to remove in order to provide the most predictive model of curvature. This process is referred to as optimization. The endpoints of optimization are the maximum coefficient of regression with the minimum VIF. Because the coefficient of regression should decrease for each curvature measure removed, there is no single answer to this problem. In order to determine the optimal subset of curvatures all possible subset combinations for a given subset size are examined through the process of an exhaustive search. The determination of the coefficient of regression for 2m subsets where m is the number of independent variables is called an exhaustive search. Multiple linear regression was performed on all 22 3-D computer models using all combinations of curvatures for a platelet radius of 3.0 mm. The results were sorted first by subset size (number of curvature measures used in the regression analysis) and then by the coefficient of regression. The curvatures that constitute the subset are those that produce the maximum coefficient of regression for a given subset size. Once the optimal subset is determined, the

VIF for each curvature in the subset is computed. The coefficient of regression, F statistic and VIF for each subset is displayed in Table II. Empty entries indicate that the curvature measure was not used in the subset.

As noted in Table II, there is no change in regression coefficient for subset sizes of eleven to sixteen subsets. Furthermore, different curvature measures can be included or excluded from an optimum subset, but RINCIR, NPER, and TILT are present as an optimal curvature measure for every subset from a size of three to sixteen curvature sets.

In an other embodiment of the invention herein, optimization of subsets can be performed by the process of backward elimination, in which each independent variable is tested to determine the variable that can be eliminated while still maintaining the largest coefficient of regression. Table II indicates that curvatures H, DK1K2V, RINCIR, NPER, and TILT are ideal for identifying aneurysms when only five curvature measures are to be used. All five of these curvature measures were computed at each vertex of all

22 3-D computer models using a platelet size of 3.0 mm.

TABLE II (continued).

At this point in the described method, the optimal curvature measures for the detection of aneurysms have been found through the optimization steps described above. The five curvature measures determined are not only optimal for the detection of aneurysms, but are also non-collinear. Therefore, each measure provides some useful information that can be used to recognize an aneurysm based on its shape. By performing multiple linear regression on these 5 curvature measures for all vertices of all 22 3D computer models in this embodiment of the invention, it is possible to determine the coefficients, j3m, of the equation

y,=0.017035 - 0.005998 H,+ 0.032173 DK1K2V, + 0.04558 RINCIR, - 0.009688 NPER, + 0.012272 TILT,. Equation 20

The above equation is optimal for detecting aneurysms using the classical and new curvature measures described. The above equation can then be applied to any future 3-D computer models of aneurysms in the following way. A 3-D, triangulated computer model of the cerebral vasculature is created by any technique (e.g., marching cubes). A bivariate quadratic polynomial is fit to the collection of vertices, i.e., the platelet, immediately surrounding a central vertex as shown in Figure 2. The subset of 5 curvature measures are calculated for each 3D surface vertex by determining their values on the bivariate quadratic polynomial. Equation 20 is then used to determine y for each vertex using the 5 individual curvature measures calculated for that vertex. A table of color values (a color table) is then created to represent the expected range of y_; and the surface vertices or surface triangles are colored according to the color table. The color table can be interactively adjustable by the user or it can be fixed. Regardless of the technique used, the intention of the color table is to provide a simple means to color the 3-D computer model such that an observer's attention is directed to an area that has curvature properties consistent with an aneurysm. The 3-D computer models can be displayed in any number of different formats including 2D black-and-white images, 2D color images, stereographic images, movie format or in 3-D computer model format. The results of the above analysis can also be used to color or enhance the original 2D cross-sectional images to highlight aneurysms on the original 2-D MRA images. While the invention has been described in terms of a single preferred embodiment, those skilled in the art will recognize that the invention can be practiced with modification within the spirit and scope of the appended claims.

Claims

Having thus described our invention, what we claim as new and desire to secure by Letters Patent is as follows:

L A method of evaluating determining three dimensional structures comprising the steps of:

a) obtaining a computerized three dimensional representation of a structure or structures; b) identifying a first set of regions on the three dimensional structure or structures and assigning a numeric value to said structure or structures; c) identifying a second set of regions and assigning a numerical value to said regions; d) determining values for a plurality of curvature measures for each vertex on a surface of the structure or structures; e) performing multiple linear regression analysis on said values determined in said determining step to obtain a coefficient of regression for all curvatures for all vertices; f) determining the variance inflation factor for each of said curvature measures; g) if all variance inflation factors are less than 10, go to step 1; h) if any variance inflation factor is greater than 10, sequentially reduce the subset of curvature measures used in multiple linear regression by 1 ; i) performing multiple linear regression on all combinations of curvature measures possible for each subset; j) selecting the subset of curvature measures that yields the largest coefficient of regression; k) performing multiple linear regression analysis on said values determined in d) to obtain a coefficient of regression for said curvature subset; and 1) inserting the partial coefficients of linear regression into the linear equation generated by multiple linear regression for said curvature subset.

2. The method of claim 1 wherein said plurality of curvature measures is selected from the group consisting of: the principal curvature kl ; the principal curvature k2; the mean curvature; the Gaussian curvature; a ratio of the minimum of the absolute value of the two principal curvatures to the maximum of the absolute values of the two principal curvatures; the difference between the first principal curvature value and the second principal curvature value; the average of the dot product between the first principal curvature vector of a surface vertex and the first principal curvature vector of all vertices on the 3-D triangulated surface immediately connected to this vertex; the average of the dot product between the second principal curvature vector of a surface vertex and the second principal curvature vector of all vertices on the 3D triangulated surface immediately connected to this vertex; the difference between the average of the dot product between the second principal curvature vector of a surface vertex and the second principal curvature vector of all vertices on the 3D triangulated surface immediately connected to this vertex and the average of the dot product between the first principal curvature vector of a surface vertex and the first principal curvature vector of all vertices on the 3D triangulated surface immediately connected to this vertex; the radius or diameter of a circle inscribed into the normal triangle such that the three sides of the normal triangle are tangent to the inscribed circle; the radius or diameter of a circle circumscribed around the normal triangle such that each of the three vertices of the normal triangle intersects the perimeter of the circumscribed circle; the ratio of the radius or diameter of a circle inscribed into a normal triangle such that the three sides of the normal triangle are tangent to the inscribed circle and the radius or diameter of a circle circumscribed around a normal triangle such that each of the three vertices of the normal triangle intersects the perimeter of the circumscribed circle; area of the normal triangle; perimeter of the normal triangle the ratio of the area over perimeter squared of the normal triangle; and the dot product of the surface triangle normal vector and the normal triangle normal vector.

3. The method of claim 1 wherein said three dimensional representation is obtained from an image of said structure.

4. The method of claim 2 wherein said image is constructed from multiple two-dimensional scans.

5. The method of claim 1 wherein said three dimensional structure is a tumor, aneurysm, or polyp, and said predetermined shape is spherical.

6. A method of evaluating three dimensional renderings of blood vessels for detecting the presence of aneurysms comprising the steps of: a)obtaining a computerized three dimensional representation of said blood vessels;; b) identifying a first set of regions on the three dimensional structure or structures and assigning a numeric value to said blood vessels; c) identifying a second set of regions and assigning a numerical value to said regions; d) determining values for a plurality of curvature measures for each vertex on a surface of said blood vessels; e) performing multiple linear regression analysis on said values determined in said determining step to obtain a coefficient of regression for all curvatures for all vertices; f) determining the variance inflation factor for each of said curvature measures; g) if all variance inflation factors are less than 10, go to step 1; h) if any variance inflation factor is greater than 10, sequentially reduce the subset of curvature measures used in multiple linear regression by 1 i) performing multiple linear regression on all combinations of curvature measures possible for each subset; j)selecting the subset of curvature measures that yields the largest coefficient of regression; k) performing multiple linear regression analysis on said values determined in d) to obtain a coefficient of regression for said curvature subset; and 1) inserting the partial coefficients of linear regression into the linear equation generated by multiple linear regression for said curvature subset.

7. The method of claim 1 or claim 6 further comprising the steps of; determining a single scalar value for each vertex on a structure or structures from the linear equation determined; creating a color table of all scalar values for each said vertex; assigning a color to each vertex; rendering the 3-D structure or structures with said color values displayed on the surface of the structure or structures.