Wire-mesh generation from image data
FIELD OF THE INVENTION 
     This invention generally relates to improvements in a computer graphic systems and more particularly to a technique that enhances a representation of image data in a graphic display system. 
     DESCRIPTION OF THE INVENTION 
     The invention converts three-dimensional image volumetric raster data set into graphic vector sets of selectable surfaces Raster data sets (images) consist of a three-dimensional array of numbers. For example, numbers indicative of intensity at each x, y coordinate on the display. Graphic data set (vectors) consist of a list of point locations and command lists indicating how the points are connected. 
     Three-dimensional raster to vector conversion is important to utilize the large variety of engineering analysis tools designed for wire mesh or engineering models, and to exploit the high speed graphic processors in presenting three-dimensional information visually. This facilitates utilization of engineering models and interactive engineering tools for visualization. 
     Current conversion tools exist for converting three-dimensional volumetric data into three-dimensional visual renderings based on perspective viewing, optical parameters of the contents, textual properties, and lighting models. These methods are directed towards perceptual quality (photo-realistic) renderings of the three-dimensional data. These techniques do not lend themselves to engineering modeling and analysis In addition, these techniques tend to be computationally intensive making interactive processing virtually impossible. 
     A classic method for generating three-dimensional vector versions of raster data is based on contouring individual planes of the data. Ad-hoc techniques are then applied to connect the contours from adjacent planes. 
     In the Computer Aided Design/Computer Aided Manufacturing (CAD/CAM) arena, methods exist for generating surface meshes for ideal geometric objects such as circles, triangles, squares, etc. and to some extent, three-dimensional objects such as spheres, cubes, and tetrahedrons. These methods do not lend themselves to arbitrary, real world objects that arise when dealing with sensory data. 
     Present image processing techniques for display on a graphic system are unable to match the original image with the precision that the invention provides Examples of various methods for converting image data into graphic representations are found in U.S. Pat. Nos. 4,879,668; 4,831,528; 4,710,876; and 4,719,585. 
     This invention is a method and apparatus to obtain well defined surface descriptions from a three-dimensional image or field data. A surface of points is modified to conform to the surface by continually refining the points until they conform precisely to the surface to be matched. The resulting set of vectors provide remarkable improvements in representing the surface to be modelled over existing methods. 
     For example, in medical imaging, it is often necessary to create a wire-frame model for three-dimensional modeling of the surface of a person&#39;s anatomy, say a person&#39;s face. For reconstructive surgery, it is imperative that the model is precise. Employing the invention, a Doctor can obtain a wire mesh that most closely approximates a piece of cellophane stretched to precisely match the patient&#39;s facial features 
     By creating a highly accurate, wire-mesh representation, current medical imaging applications, Magnetic Resonant Imaging (MRI) techniques, CAT Scans and existing CAD/CAM analysis routines can be more effectively exploited. 
     SUMMARY OF THE INVENTION 
     An integrated graphic and image system for generating a wire-mesh representation of graphic information representative of a three dimensional surface in image data is disclosed. The system includes a disk drive for storing graphic information representative of a three dimensional surface. The system also has a memory and central processing unit (CPU) for generating data representative of an initial wire-mesh comprising a set of points in three dimensional space and a topology. The CPU further refines the initial wire-mesh representation by modifying the initial wire-mesh to intersect the graphic information representative of the original three dimensional surface at specific points. The CPU generates a new wire-mesh by clamping the intersecting points to the image information representative of a three dimensional surface to further refine the wire mesh representation. The CPU repeats the refining steps until every point of the initial wire-mesh is clamped to form a wire-mesh representation of the three dimensional surface. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a block diagram of a computer graphic system in accordance with the present invention; 
     FIG. 2 is a flow chart in accordance with the present invention; and 
     FIG. 3 is a time diagram in accordance with the present invention. 
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     Referring first to FIG. 1 of the drawings, there is shown a block diagram of a computer system employing the invention. The logic represented in the flowcharts (FIG. 2) resides in the memory 30 of the graphic system 10. A CPU 20 communicates with the memory 30 and executes the instructions attached below in the listing of the program in accordance with the present invention. Various other units can be connected to the CPU 20 in addition to the graphic display 50 and the disk 40. 
     In practice, the program is loaded from the disk 40 and reads the image data and generates the initial mesh as depicted in function block 100 of FIG. 2. The initial mesh is generated as a surface as shown at time zero in FIG. 3 at label 200 on the graphic display 50. The initial mesh is composed of mobile masses 250 connected by connecting springs 270. The mobile masses and connecting springs are represented in the memory of the graphic system by mathematical equations predetermined by the user based on the characteristics of the surface to be represented. 
     An embodiment of the invention currently resides on an IBM 3090 computer. The display is an IBM 6090 (or 5080) graphic workstation. The computer system comprises a keyboard, display, memory, disk and S/370 processor. The graphic workstation comprises a local memory, display and various input devices. The input devices include a keyboard, tracking cursor with tablet, dials and a multi-function keypad. IBM&#39;s graPHIGS software is used to create the images and refine the surfaces in accordance with the present invention and the program listing which appears below. 
     Various examples of effective equations are discussed below as examples utilized currently. However, the invention is not dependent on the use of any particular equation. After initiation, the initial mesh is moved into intersection with the surface to be represented, depicted as a hump at the base of FIG. 3 at label 240. The mobile masses 250 become fixed masses 260 at the points of intersection with the surface as depicted at time=two 220 of FIG. 3. The actual positions of the masses and their points of intersection are updated dynamically on the display as the initial mesh is moved into place as shown in function block 120 of FIG. 2. However, the mobile mass 250 is not transformed into a fixed mass 260 until a trajectory segment intersects the surface to be matched. A trajectory segment refers to the segment between two mobile masses. 
     The refinement of the initial mesh by moving it into more and more of a one-to-one correspondence with the surface to be matched continues until all of the mobile masses have been transformed into fixed masses as represented in decision block 140. When the final mesh is complete, it is displayed on the display 50 and written to disk 40 as shown in function block 150. 
     The apparatus and method discussed above processes three-dimensional volumetric raster data and produces a graphic vector data set. This output graphic vector data set is a set of three-dimensional points and associated connectivity between those points that produce a wire mesh rendering of a selected surface in the three-dimensional raster data. The volumetric data sets can be arbitrarily complex--no simple geometric object decomposition is assumed. 
     The three-dimensional points and the associated connectivity is referred to as a geodesic. 
     
         G={P.sup.i,T.sup.i }i=1,N 
    
     P i  =(x i ,y i ,z i )=vector data set of node points on surface 
     x i  =x coordinate of point i 
     y i  =y coordinate of point i 
     z i  =z coordinate of point i 
     T i  =topology list=list of neighbor nodes for node i 
     Given a three-dimensional multicolored image, we seek to find a geodesic representation of surfaces in the image. Let I xyz  =source image of colors located in three-dimensional space. 
     In this case, a surface in the image means an exterior boundary for which the color (or intensity) is a constant: 
     
         G={P.sup.i,T.sup.i } 
    
     I p  =I O  =a given threshold constant 
     Note that the geodesic format can easily be converted to other wire mesh specifications such as facets models or simple move/draw command lists. The derivation given here is based on the above format; conversions to other forms is self evident. The simplest wire mesh rendering is to draw line segments between every point and all of its neighbors: 
      i draw line from P i  to P j  where jεT  i   
     The algorithm emulates a &#34;shrink wrap&#34; process by simulating the contraction of an elastic geodesic. The simulation is performed by mathematically modelling a hypothetical elastic membrane undergoing a collapse due to internal suction or from surface tension. Variations of the basic method allow for different dynamic models for the collapse, different physical parameters for the elastic nature, and different initial conditions for the starting mesh. 
     SHRINK WRAP 
     The method and apparatus contract a given geodesic to the image surface. The contraction is performed via a time simulation. Assume we are given an initial geodesic consisting of node locations and an associated topology: 
     
         G(t=0)={P.sup.i,T.sup.i } 
    
     where t=time 
     The initial selection of the geodesic is a system input. For clarity, assume for now a longitude/latitude type grid with equal angular spacing. The points (P) are then the three-dimensional locations of the grid intersections and the topology (T) for each point refers to each of the four neighbors (above, below, left, and right). Without loss of generality, assume that the image lies entirely within starting geodesic set (if not, simply magnify the geodesic by the size of the image). 
     The output geodesic is computed by assigning physical parameters to each of the nodes and topology links and then allowing the resulting physical system to transition to equilibrium. This physical model is what is commonly known as a lumped parameter system. Each node is assigned a mass and a model of the forces is established for each of the topology links. 
     Topology link forces can be extremely simple (e.g. free fall), linear (e.g. Hooke&#39;s law for springs) or non-linear (e.g. terminal velocity) all of which are described in more detail below. Whatever physical model is selected, the initial geodesic is collapsed onto the surface by simulating the Newtonian physical model based on position P(t), velocity P(t), and acceleration P(t). (Super &#34;dot&#34; denoting time derivative.) ##EQU1## 
     The time simulation is performed by numerical integration; the simplest means being to use Eulers method although other techniques are equally applicable (such as Runge-Kutta). Eulers method for a given δt becomes: ##EQU2## 
     At each time update, a check is made to determine if the trajectory from P i  (t-1) to P i  (t) intersects the surface of interest. This test is performed by testing the intensities at all unique image locations along the trajectory and determining if and when the intensity profile crosses a given threshold. Resampling the image data at trajectory positions is performed using any of the standard methods (e.g. nearest neighbor, tri-linear interpolation, etc.). If the threshold is crossed (i.e. if the trajectory intersects the surface) then the associated node position is set to the surface location. In addition, the mass associated with that node is set to infinity. This has the effect of fixing that particular point to the surface and no further motion of that node will occur. 
     if [P i  (t-1), P i  (t)] intersects surface 
     P i  (t)=surface point 
     M i  =∞ 
     (In actual implementation, rather than setting the mass to infinity, the node is simply tagged as being fixed and the position and velocity updates are no longer performed.) 
     The time integration proceeds iteratively until all motion has ceased or until the total kinetic energy has decreased below a threshold. The desired geodesic is given by the steady state (static equilibrium) result: ##EQU3## 
     Static equilibrium can occur when all points have in fact collapsed onto the surface or when the geodesic has been stretched across open regions. Points that are not on the surface can be pruned out from the geodesic simply by deleting the point and all references to that point in the topology list. 
     Static equilibrium will occur provided there are dissipative forces and that the external forces are designed to drive the system to a fixed point. 
     SPECIFIC MODELS 
     For completeness, specific physical models are now presented. The pros and cons of the models are subjective and case dependent. In fact, a whole plethora of physical models exist to model different types of elastic material in a variety of environments. The method presented proceeds using numerical integration of the system state which easily encompasses any physical model of the general form: 
     
         F.sup.i (t)=f(P.sup.i (t-k), P.sup.i (t-k), P.sup.i (t-k), t) 0&lt;k≦t 
    
     SPECIFIC MODEL: FREE FALL (FF) 
     In this model, FF, topology link forces are non-existent. An internal suction (vacuum) force is applied. This force is such to draw all node points towards a common origin. In addition, a dissipative friction force is applied to force the system to converge. In this case, the shrink wrap method reduces to projecting rays from each of the initial geodesic locations into the image and determining the resulting position. For given constants F 0  and b i , 
     
         F.sup.i (t)=-F.sub.0 P.sup.i (t-1)-b.sup.i P.sup.i (t-1) 
    
     SPECIFIC MODEL: LINEAR LUMPED PARAMETER (LLP) 
     The LLP model assumes an elastic model based on idealized linear springs (Hooke&#39;s law) and a dissipative energy proportional to the velocity. Each link is assigned a spring constant (K ij ) and each point is assigned a frictional constant (b i ). ##EQU4## 
     MODEL VARIATION: TERMINAL VELOCITY 
     Both of the above models can be varied by altering the dynamic model. One such variation is to hard limit the velocity giving FF-TV and LLP-TV models. Thus, the acceleration and position updates are as before but the velocity equation becomes: ##EQU5## 
     OTHER MODELS 
     Non-linear extensions to Hooke&#39;s law and other non-linear versions of frictional dissipative forces are easily incorporated into the scheme. The above models described above disclose techniques that varied the physical model and parameterization resulting in different state trajectories on P. Another variation that the proposed technique allows is to dynamically alter the topology. For example, elastic breaking can be simulated by deleting neighbors of a node. Thus, at each time step, for a given break threshold d break , ##EQU6## 
     INITIAL GEODESIC 
     The initial geodesic is a predefined surface initially set at system definition. The accuracy of the resulting geodesic is directly correlated to the surface chosen for the starting geodesic. Factors that yield favorable results on actual objects are presented. 
     1) In general, initial geodesics should correspond to the gross approximate shape of the object. 
     2) A spherical longitude/latitude grid with uniform angular increments generates excellent representations for medical data of a human skull. 
     3) A cylindrical grid with uniform angular increments and uniform spatial increments along the z axis was found to generate excellent representations for certain industrial parts. 
     4) Non-uniform Rational B-Splines (NURBS) are a well known surface representation technique that can be exploited by selecting a NURB surface as the initial geodesic of the designers choice. 
     Non-Uniform Rational B-Splines (NURBS) are used in one embodiment to describe all large complex surfaces mathematically. As described in IBM&#39;s AIX Personal graPHIGS Programming Interface/6000 Version 2, copyright Jan. 29, 1990, NURBS are supported by IBM in their graphic programming language. NURBS and the subsequent productivity boost incorporated by their use is described in the October 1989 Computer Graphics REVIEW on pages 22 to 26. NURBS support complex mathematical definition of surfaces in a single data structure. 
     The NURB processing manipulates a set of three-dimensional control points arranged on an ordered two-dimensional grid. In accordance with the subject invention, the ordered two-dimensional grid specifies the topology of the geodesic while the three-dimensional points are the nodes. The resulting geodesic is a NURB surface whereby all of the control points (nodes) lie on the actual surface. Further, the topology remains unchanged. Thus, the resulting geodesic is a NURB representation of the image object that corresponds very closely to the original image. 
     A computer program listing in accordance with the invention is found in the Appendix. ##SPC1##