WO1994024640A1 - System and method for surface rendering of internal structures within the interior of a solid object - Google Patents

System and method for surface rendering of internal structures within the interior of a solid object Download PDF

Info

Publication number
WO1994024640A1
WO1994024640A1 PCT/US1994/004084 US9404084W WO9424640A1 WO 1994024640 A1 WO1994024640 A1 WO 1994024640A1 US 9404084 W US9404084 W US 9404084W WO 9424640 A1 WO9424640 A1 WO 9424640A1
Authority
WO
WIPO (PCT)
Prior art keywords
voxel
vertex
surface structure
face
intersect
Prior art date
Application number
PCT/US1994/004084
Other languages
French (fr)
Inventor
Daniel Karron
James Cox
Bhubaneswar Mishra
Original Assignee
Daniel Karron
James Cox
Bhubaneswar Mishra
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Daniel Karron, James Cox, Bhubaneswar Mishra filed Critical Daniel Karron
Publication of WO1994024640A1 publication Critical patent/WO1994024640A1/en

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20Finite element generation, e.g. wire-frame surface description, tesselation

Definitions

  • This invention generally relates to a system and method for surface rendering and more particularly to a system and method for the display and analysis of arbitrary three dimensional surface structures.
  • a typical application is in medicine where the solid object can be a human body or a part of it and the objective is to detect and model certain internal structures for diagnostic and surgical planning purposes.
  • the methods to which the present invention relates are, however, equally applicable to other natural or artificial bodies and fields of study.
  • Examples of the physical properties measured at regular 3D positions include the coefficient of x-ray absorption in the CAT case or, for NMR imaging, the spin-spin or the spin-lattice relaxation parameters.
  • the measured physical values reflect variations in the composition, density or structural characteristics of the underlying physical objects, providing knowledge about their internal structure.
  • the 3D data array typically corresponds to a plurality of measurements of the given physical property, taken at regular intervals along an appropriately selected system of 3D (X, Y, Z) physical coordinates.
  • Each measurement V ⁇ of the physical property is thus associated with a specific array index position in a cubic lattice within the volume in study.
  • a cubically adjacent set of eight such positions defines an unit cubic volume of the body under study and is known in the art as a voxel.
  • bones, brain or any other tissues can be characterized by the known range of their density values which is compared to the measured and stored 3D array values in order to separate out a surface model of their structure.
  • the displayed visual images can next be used to assist the medical practitioner in providing a diagnosis or in the planning of complex surgical procedures.
  • this surface is typically illuminated and shaded so as to give the human eye the impression of a natural shape and spatial disposition of the surface when it is displayed.
  • the angular direction of a vector orthogonal to the surface at each point on the surface is compared to the viewing angle of the observer.
  • the intensity of shading can be adjusted to be proportional to the difference between these angles.
  • the surface segment intersecting a voxel is approximated by one of a limited number of standardized plane polygonal surfaces intersecting the voxel.
  • One particular standardized surface tiling is selected from a library of available surfaces by a vector which represents the binary differences between the threshold value for the given surface of interest and the eight voxel vertex values.
  • the surface-to-voxel intersection coordinates, as well as the vector normal to the surface for each such standardized polygonal surface tiling can be obtained by table look-up techniques.
  • the final surface is assembled as a mosaic using all the standardized triangular tiles.
  • the values at the vertices of the voxel are used to interpolate in three dimensions values at regularly positioned intra-voxel sub-grid locations.
  • the interpolated sub- grid values can then be used to locate the surface position and calculate the normal vector more accurately.
  • U.S. Pat. No. 4,751,543 to Lorensen discloses a technique for labeling surfaces with similar property values and using adjacency criteria in respect to a "seed location" in the particular structure of interest to segregate the desired surface from all of the labeled surfaces.
  • U.S. Pat. No. 4,791,567 to Cline et al. discloses another technique of segregating similar internal structures by determining connectivity from adjacency information.
  • U.S. Pat. No. 4,821,213 to Cline et al. discloses yet another technique for differentiating internal structures in which a linear path is made through the data array to locate and label all different structures along the scan line by counting structure interfaces.
  • voxels within the solid object under study are typically traversed sequentially and progressively, so that every voxel within the data space is checked for a crossing by the surface of interest, and is subsequently assigned to belong either to the outside or to the inside of the surface. Because of the large volumes of data to be processed, the methods are computationally intensive and relatively slow. An improvement in this respect is described in U.S. Pat. No. 5,166,876 to Cline, issued in November 1992. The disclosed method significantly speeds up the rendering of a surface within a boundary by only processing voxels on the surface. This approach saves significant computation time but still uses the same case table and the limited set of approximating surfaces which were described in the original marching cubes patent.
  • the marching cubes surface rendering method is widely accepted but has certain deficiencies. These deficiencies are generally related to the fact that the fifteen canonical cases used in the surface segments' approximation do not cover the entire spectrum of possible tiling combinations. This results in occasional approximation errors, occurring at lower levels of resolution, specifically if a voxel face is intersected by the surface twice resulting in four crossings ("hits") with the voxel edge. The errors manifest themselves in incorrect surface rendering, visually appearing as unexpected spots, tears and missing tile elements. The problems are aggravated if the collection data is noisy and/or with low contrast. A number of corrections to the marching cube algorithm have been offered in the past.
  • the display of 3D graphic images has been principally driven by the goals and directions of computer aided design (CAD) and computer aided manufacturing (CAM) systems.
  • CAD computer aided design
  • CAM computer aided manufacturing
  • Systems accepting data in several basic formats have been developed for displaying solid bodies and manipulating images in various ways to create models of parts being manufactured, or internal organs being monitored.
  • the present systems for surface rendering and display are, however, still incapable of meeting the requirements for a fast and truly interactive display of continuously monitored 3D data.
  • the complexity of the task can be better appreciated by noting that a typical surface rendering operation involves processing on the order of 512 x 512 x 512 pixels. In the presence of noise the task becomes even more complicated because of tbe necessity for additional noise filtering, smoothing and data interpolation, which takes processing time and often results in the loss of important structural details.
  • a further object of the present invention is to provide a system for analysis of internal structures within a 3D solid body.
  • a novel method for rendering three- dimensional surface structures comprises the following steps.
  • Signal patterns representing the value of at least one physical property associated with regularly spaced grid locations in a three dimensional body are stored in a memory means.
  • the regularly spaced grid locations define voxel elements within the body, such that each voxel is associated with eight vertex values and six voxel faces.
  • the signal patterns that represent values lying within a specified range associated with the surface structure of interest are next selected for further examination.
  • An initial grid location which corresponds to a position within or on the surface structure is first selected and the eight vertex values associated with a voxel at the initial grid location are selectively retrieved.
  • a threshold value is supplied to determine which voxel vertex values are inside and which voxel vertex values are outside of the surface structure. For each voxel face having at least one voxel vertex inside and at least one voxel vertex outside of the surface structure, there must be at least one pair of intersects between the edges of the voxel face and the surface structure. The coordinate values of these intersects between the voxel face and the surface structure are then computed. Preferably, the computation of these intersect coordinate values is done by a linear interpolation of the voxel vertex values. Other methods of computation can similarly be used.
  • a triangular approximation to the surface structure is then formed using the pair of intersects between the edges of the voxel face and the surface structure as two vertices of the triangle with a third vertex in the interior of the voxel.
  • the step of forming triangles is repeated for each pair of voxel face intersects on the six faces of the voxel.
  • the physical location of the third vertex of each triangle is identical and is computed as the centroid of all interpolated intersect coordinate values for the voxel.
  • the method further comprises the steps of identifying all adjacent voxels sharing a face intersect with the surface structure.
  • the steps of retrieving the eight vertex values, determining the inside and outside vertex values, computing the coordinates of the voxel face intersects and forming the triangles is repeatedly applied until all face-adjacent voxels within the three-dimensional body having a face intersect with the surface structure are processed.
  • the step of computing the intersect coordinates only values that have not previously been computed are processed.
  • the proposed method of surface rendering further comprises the step of displaying a two-dimensional image of the approximated surface structure.
  • the surface structure is displayed by utilizing only the voxels identified by the processor as being crossed by the surface structure.
  • a method for generating of an internal surface structure within a three-dimensional body comprises the following steps.
  • Three-dimensional signal patterns representing the value of at least one physical property associated with regularly spaced grid locations in a three dimensional body are stored in memory means.
  • the regularly spaced grid locations define voxel elements within the body, such that each voxel is associated with vertex locations and voxel faces, each two intersecting voxel faces defining a voxel edge.
  • Each voxel edge passes through two voxel vertices.
  • the signal patterns which represent values lying within a specified range associated with the surface structure of interest are next selected. Two voxel vertices of a voxel edge for one of said stored voxels are then retrieved; and it is determined which voxel vertex of the selected voxel edge is inside and which voxel vertex of the selected edge is outside of the surface structure.
  • a coordinate value of a first hit location where the surface structure crosses the voxel edge is next computed.
  • the computation of the intersect coordinate values is done by a linear interpolation of the voxel vertex values but other computation methods can be used as well.
  • each of the two intersecting voxel faces that define the voxel edge where a hit was identified is examined for a second hit location which defines an intersect between the surface structure and the voxel face.
  • a triangle with a vertex in the interior of the voxel is formed to obtain an approximation to the surface structure.
  • the location of the third vertex of each approximating triangle is computed as the centroid of the voxel face intersects such that the voxel face intersects form a closed contour around the inside voxel vertex.
  • the method further comprises the steps of identifying all adjacent voxels that share a face intersect with the surface structure.
  • the steps of selectively retrieving two voxel vertices of a voxel edge, determining the inside and outside vertex values, computing coordinates of the first hit locations, examining the voxel faces defining the voxel edges, and forming the approximation triangles are repeatedly applied until all face-adjacent voxels within the three-dimensional body having a face intersect with the surface structure are processed.
  • the step of computing the intersect coordinates only values that have not previously been computed are processed.
  • the method further comprises the step of displaying a two- dimensional image of the approximated surface structure.
  • the surface structure is displayed by utilizing only the voxels identified by the processor as being crossed by the surface structure.
  • the proposed method for 3D surface rendering is implemented on a system for analysis and display of three-dimensional surface structures which comprises the following elements.
  • storage means are provided for storing 3D signal patterns representing the value of at least one physical property associated with a 3D body at regularly spaced grid locations within the body.
  • the signal patterns are provided by a three-dimensional raw data source, such as MRI, CAT, x-ray, ultrasound or others.
  • a set of 3D-adjacent grid locations define a volume element referred to as voxel having associated with it voxel edges and voxel faces.
  • the system also comprises means for defining a plane of pixel positions having preselected values and corresponding to one plane of the spaced grid locations. Also provided is means for determining in the defined pixel plane, a pixel associated with a voxel that is intersected by a surface structure of interest within said body. The system further comprises processing means for computing coordinate values of the intersection between the surface structure of interest and the intersected voxel edges and for computing the coordinate values of at least one articulating location positioned in the interior of the intersected voxel. Second storage means for storing the computed coordinate values along with connectivity information about joining the computed coordinate values into a set of triangles is also provided. These triangles approximate the intersection between the surface structure of interest and the voxel.
  • Each triangle has one edge on a voxel face and one vertex at an articulating location in the interior of the intersected voxel.
  • the processor means further comprises generating means for generating normal vectors associated with the computed coordinate values.
  • the computed normal vectors are similarly stored in the second storage means.
  • the processing means of the system further comprises means for identifying which of the voxels sharing a voxel face with the intersected voxel are intersected by the surface structure of interest.
  • the processing means further comprises means for generating, for each voxel that is intersected, a plurality of index values representing information about the shared voxel face intersects. These index values are also stored in the second storage means and are used to indicate to the processor which voxels and voxel faces have been processed.
  • the processing means of the system comprises means for simultaneously comparing all vertex values of an intersected voxel to a threshold value. This comparison is done independently for each voxel face thus enabling a parallelizable, high speed implementation of the proposed system.
  • the system of the present invention further comprises display processor means which receives data from the second storage means and converts it into a display format.
  • the display format data is supplied to means for interactively displaying surface structures of interest.
  • the planar arrangement of the regularly spaced grid-locations is other than a rectangular array.
  • the planar arrangement of the grid locations may have a hexagonal, or other pattern suitable for processing shape.
  • the system and method for surface rendering and display of the present invention have a number of practical applications.
  • the internal data structure representation allows the user to simultaneously process and display two or more isosurfaces.
  • the system provides a mechanism to count 3D disjoint structures within the interior of a body under study by enumerating the independent structures being processed.
  • the system and method of the present invention can in addition be used to compute and display cross- sections of disjoint surface structures.
  • Such surface structures need not only be internal structures of the object under study. While the description of the proposed method is presented in terms of the measured values of a given physical property of the object, other types of surfaces can similarly be computed and displayed. For instance, this can include isogradient and isocurvature surfaces giving information about the locations within the body having identical rate of change in the density values or curvature values, respectively.
  • the cross-section of such surfaces with the structure surface can be used to automatically compute cardinal landmarks on the object surface.
  • the system of the present invention can be used to compute global structure parameters such as the total surface area, and the structure volume. Information necessary to compute these quantities, including the coordinates of a set of physical points which belong to the surface is stored and available for each processed structure. The resulting parameter estimates can be used to define characteristics of the object under study, such as volume occupancy, average structure size and others, which characteristics have a wide range of applications, especially in 3D pattern classification.
  • the system of the present invention can be used to observe time changes in the interior structures of a given solid object by knitting together "snap-shots" of the surface of the structure and displaying these snap-shots in a sequence. An exemplary application of such system use is the monitoring of a patient's heart beat pattern.
  • Fig. 1 is a detail of surface tiling errors in a prior art surface rendering method.
  • Fig. 2 is a block-diagram of a system for carrying out the method of the present invention.
  • Fig. 3 is an isocontour plot of an MRI axial slice of a brain used for plane selection preview.
  • Fig. 4 is a three dimensional voxel with elements defined.
  • Figs. 5(A) and (B) are perspective views of a voxel neighborhood and voxel adjacencies for three dimensional voxel elements, respectively.
  • Fig. 6 shows intersection of a voxel by surface structures.
  • Fig. 7 illustrates voxel face hits assumptions.
  • Fig. 8 shows ambiguity in voxel face contour lines.
  • Fig. 9 is a flow-chart of the SpiderWeb A method.
  • Fig. 10 illustrates spinning approximation triangles using an articulation point.
  • Fig. 11 illustrates spinning approximation triangles for a voxel face contour ambiguity using SpiderWeb A.
  • Fig. 12 is a flow-chart of the SpiderWeb B method.
  • Fig. 13 shows a two-dimensional adaptive sampling grid.
  • Fig. 14 shows triangle vertices drawn in a T-mesh order.
  • Fig. 15 is an approximating surface computed with full precision.
  • Fig. 16 shows the approximating surface of Fig. 16 using midway hit location speed processing.
  • Fig. 17 shows triangle face and vertex normal vectors.
  • Figs. 18 and 19 illustrate the medical imaging surfaces rendered in accordance with the method of the present invention.
  • FIG. 1 a surface detail is shown which illustrates tiling errors in a prior art surface rendering method.
  • the errors are observed as bright or very dark spots on the surface, indicating a tear in the represented surface or a missing approximation tile.
  • the present invention is designed to remove such errors by presenting a topologically consistent method of building a triangularization grid which represents an arbitrary three-dimensional structure.
  • the invention is practiced on the raw data obtained from non-intrusive measurements of at least one physical property within a planar slice through a body under study.
  • a typical source of the raw data is an MRI or CT system, an x-ray, ultrasound or another type of a scanner, capable of converting measurements of some physical property of the body, at regular grid-like locations within the planar slice, into an array of discrete values.
  • Other non-intrusive data collection mechanisms such as positron emission tomography (PET) , emission computer tomography (ECT) and multimodality imaging (MMI) systems can also be used as a data source.
  • PET positron emission tomography
  • ECT emission computer tomography
  • MMI multimodality imaging
  • the operation of the data source may be carried out for a sequence of slices of the solid object, thus providing a three-dimensional array of values representative of the spatial structure of the object.
  • the raw data is stored in a first storage means, preferably indexed by a slice number.
  • the storage means preferably has a fast access time and may comprise magnetic storage, random access memory or other known in the art storage media.
  • the storage means also comprises means for converting the raw data into an application specific internal data structure, in order to take advantage of the specific hardware configuration used.
  • the internal data structure provides the capability of individually addressing and retrieving individual pixel elements as well as sets of spatially adjacent pixel configurations, referred to as voxels.
  • the first storage means also provides the capability of indexing each vertex, each voxel edge and voxel face, in order to facilitate the computation process.
  • Fig. 4 shows a voxel element consisting of a cubic arrangement of eight data pixels and presents element definitions consistently used in the following description.
  • Each pixel element is associated with a triplet of indexes (i,j,k) used to determine the location of the pixel in a three-dimensional indexing system.
  • indexes i,j,k
  • the lower left vertex of the front face is selected as a reference point.
  • the indexes of the adjacent pixels to the left, to the top and towards the back plane are shown to increase by the length s, which corresponds to the number of display pixels inscribed within two consecutive data points. It is to be understood .that Fig.
  • voxel elements such as the vertices, edges and voxel faces, illustrated in Fig.
  • a second storage means associated with the data processor comprises a memory which is allocated for each data structure, so that the processor is capable of determining whether or not a particular element of the data structure has been previously processed. For instance, each time a specific voxel edge is retrieved, a mark is left indicating for all associated data structures that this element has been processed. Thus, when the next element of a connected data structure is visited, the processor need not revisit the original voxel edge.
  • Fig. 5 (A) is a perspective view of a voxel neighborhood.
  • Fig. 5 (B) depicts possible combinations for voxel pairs, which define a vertex adjacency, an edge adjacency and a face adjacency, depending on whether only a vertex, an edge or a voxel face, respectively, are shared.
  • a voxel element is shown intersected by two separate isosurfaces. Voxel vertices, having values above an isosurface threshold are termed "inside" vertices, while voxel vertices having values less than the isosurface threshold are termed "outside" vertices.
  • a hit The intersection of an isosurface of interest with a voxel edge is referred to as a "hit”.
  • Six such hits, indicated as H0-H5 are shown in Fig. 6.
  • the intersection between the surface of interest and a voxel face is referred to as a "face intersect” and, as seen in Fig. 6, connects two voxel edge hits.
  • the first storage means comprises means for converting the raw data into a specific internal data structure, which allows all voxel elements of the data array to be consistently indexed by the data processor for identification purposes.
  • the face intersect H0-H1 on the left face of the voxel in Fig. 6 is shared with the adjacent left voxel and, in fact, occupies the same physical location.
  • the data processor need only compute the hit locations once, and use their values when processing all corresponding edge adjacent voxels, as defined in Fig. 5(b).
  • the data processor implements a topologically consistent method of connecting the isosurface hits on voxel edges into a continuous triangularization grid representing the surface of interest for analysis and display.
  • the method requires a starting point, or a "seed" pixel element on a data plane, selected by a plane selecting means.
  • the method uses only local information about the adjacent pixels, the method "spins” a web of triangles, closely approximating the surface of interest.
  • the processing sequence and visual appearance of the rendered surface resemble the construction of spider web, giving the process its "Spiderweb" name.
  • FIG. 7 depicts a voxel face with vertices V0-V3, edges E0-E3 and a clockwise scanning direction, as indicated by the arrow.
  • a hit HI is shown on voxel edge El.
  • an isosurface is assumed to have no singularities, breaks, or holes, so that given a hit on one edge, there will be at least one more hit on the same voxel face.
  • hits on a voxel face must come in pairs. Referring back to Fig. 7, this implies that there must be at least one additional hit on the voxel face. If no i such hit, or an odd number of hits is found, the processor will indicate an error.
  • Figure 8 illustrates another possible situation, when all four edges of a voxel face have isosurface hits. As shown, this is only the case when diagonally opposite voxel vertices are both "inside,” or both “outside” the surface of interest, and the two diagonal pairs are not on the same side of the surface. Lacking global information about the isosurfaces, the exact direction of surface intersects cannot be determined. This ambiguity results in four possible cases of face intersects, indicated as cases A, B, C and D in Fig. 8. Cases B and D will be re ⁇ considered in connection with the SpiderWeb A, and SpiderWeb B processing methods described next.
  • Figure 9 is a flow-chart of a first embodiment of the surface rendering method in accordance with the present invention, referred to as SpiderWeb A.
  • this method no decision about the exact direction of the face intersects in the case of four hits is made. This corresponds to case D in Fig. 8, where the four hits are connected, thus minimizing the probability that an error is made.
  • the algorithmic steps in this case are as follows.
  • step 200 of the method the pixel elements corresponding to one voxel are simultaneously retrieved from the data storage.
  • the vertex elements f the voxel are next compared in step 210 to a preselected threshold which corresponds to the surface of interest, and each vertex is assigned as either inside or outside with respect to the surface.
  • the method proceeds to retrieve a face adjacent voxel in accordance with predefined voxel sequencing rules. If a surface intersection is detected at step 220, the locations of all edge hits of the voxel, which have not been previously processed, are computed. Preferably, the position of a hit along an edge is computed by using a linear interpolation of the values at the two vertices defining the edge.
  • a voxel edge with a hit has an outside vertex at position V x0 , with a value of V 0 , and an inside vertex at position V xl with a value of V, and the threshold value is T, SO that V 0 ⁇ T ⁇ V, the location of the edge hit Ex can be calculated using the formula:
  • the next step 240 involves the calculation of an Articulating Point (AP) which is located entirely in the interior of the voxel.
  • AP Articulating Point
  • the articulating point in the SpiderWeb A method is computed as the mean (centroid) of the system of all connected hits in a voxel.
  • the articulating point is used in the next step 250 as the third vertex of a set of surface approximation triangles, each triangle being defined by the articulating point and each of the possible different pairs of the connected hits in the voxel.
  • a set of surface approximation triangles each triangle being defined by the articulating point and each of the possible different pairs of the connected hits in the voxel.
  • the articulating point will lie on the plane defined by the hits; and the approximation triangles will be co-planar.
  • An illustration of the SpiderWeb A two-face-hits approximation surface is shown in Fig.
  • step 260 of the method triangle and vertex normals for each of the appr.oximation triangles are computed. The purpose of these normals will be discussed in more detail below in connection with the description of display processor 70 of the system of Fig. 2.
  • step 270 for each intersected voxel face, the face adjacent voxel is retrieved and processed. In accordance with the present invention, the indices of all voxel faces that have already been processed are marked so that processed voxel faces need not be considered again.
  • SpiderWeb A processing method has the disadvantage of creating "tunnel” voxels and "bubbles" in the case of a voxel face ambiguity.
  • this property may actually reduce the visual impact of the noise, since its effect is localized to a particular voxel element and is not visible to the user.
  • FIG 12 is a flow-chart of the second embodiment of the surface rendering method in accordance with the present invention, referred to as SpiderWeb B.
  • SpiderWeb B a consistent decision is made about the exact direction of the face intersects in the case of four hits on a voxel face.
  • the decision corresponds to case B in Fig. 8, where the four edge hits are connected in a way so as to surround the inside voxel vertex.
  • the algorithmic steps in this case are as follows.
  • step 300 of the method a starting or seed vertex is selected.
  • the three edges associated with the seed vertex are scanned so as to retrieve the vertex elements at the opposite ends of those edges, and in the following step 320 those vertex elements are compared to the pre-selected threshold.
  • the hit locations are calculated by means of a linear interpolation, as discussed above, or estimated in any other suitable manner, as known in the art.
  • Triangles are formed in SpiderWeb B in the following four steps which are performed for each edge that is hit.
  • step 340 the two voxel faces that define the edge that has been hit are scanned for additional hits. This is done beginning with the inside vertex of the edge that has been hit and proceeding from vertex to vertex around the edges of the first voxel face and then around the edges of the second voxel face. Hit locations are calculated by linear interpolation as discussed above. Since face intersects on a voxel face occur in pairs, the scanning of each face stops when the first additional hit is located.
  • the processor then forms a triangle having one vertex at the intersection point on the edge that was originally hit, a second vertex at the location of the first additional hit and a third vertex at an articulating point with an initially undetermined physical location inside the voxel. If the foregoing procedure identifies any single hits on another face of the voxel, that face is likewise scanned to locate the second hit; and the two hits on that face and an articulating point at an undefined location within the voxel constitute another surface approximation triangle. If a pair of hits are identified on another face of the voxel, these hits and an articulating point likewise constitute another surface approximation triangle. This process .
  • step 380 For each voxel edge with a surface hit, the process repeats itself in the face adjacent voxels.
  • step 380 a set of edges from the face adjacent voxel is examined until all possible 5 intersections of the isosurface of interest are processed.
  • step 390 the processor exits to storage means 60 and the display processor 70.
  • this embodiment of the present invention does not necessarily involve o simultaneous processing of entire voxels. Instead, after the initial seed pixel is specified, the proposed method tracks along the isosurface of interest in accordance with the adopted indexing scheme, only visiting voxel edges and faces with hits. 5 Both SpiderWeb A and B processing methods reflect the uncertainty of the voxel face contour line ambiguity, illustrated in Fig. 8. In the embodiment of the present invention, using SpiderWeb A, no assumption is made about the exact direction of the 0 face intersects, resulting in a three dimensional "bubble". In the second embodiment of the present invention, using SpiderWeb B, a systematic choice to "surround inside vertex" is made based on physical considerations.
  • Figure 13 shows a two-dimensional adaptive sampling in which the upper left corner is illustratively subdivided until adequate resolution is achieved.
  • sampling elements considered thus far is only illustrative of the different sampling grids which can be used in accordance with the present invention.
  • Well known in the literature is for instance the two-dimensional hexagonal sampling grid that can reduce the number of required samples in a plane by about 13.4% (see, for example Dudgeon et al., "Multidimensional Digital
  • sampling grids can also be used, as known in the art.
  • illustrative apparatus for the practice of the invention comprises a first storage means 20, a plane selection means 30, control unit 50, a Spiderweb processor 100, a second storage means 60, a display processor 70 and a structure analysis means 80.
  • Plane selection means 30 operates under the control of an operator from control unit 50 to select data from individual data planes stored in data storage 20.
  • the operator may preview data planes corresponding to a particular physical location of interest in the body under study. Since structures having different physical properties exhibit different ranges of measured values, the operator can isolate a specific structure of interest by selecting a particular threshold or a range of threshold values. For example, the operator could isolate bones from the surrounding soft tissue in medical imaging.
  • Figure 3 illustrates an isocontour plot of an MRI axial slice of the brain used for a plane selection preview in accordance with a preferred embodiment of the present invention.
  • the operator can choose a specific plane of interest and display it on a monitor.
  • a preselect threshold value is supplied to display processor 70 and is used to display the locations (the contour lines) at which the selected data plane is crossed by a surface structure having the same physical property value.
  • the operator can then identify a particular contour line and specify the coordinates of a desired starting point ("seed") .
  • the coordinates of the seed are next supplied to processor 100 to initiate the rendering of the desired surface structure. If no contour line of the surface of interest is identified, the operator can interactively specify a different data plane or modify the threshold value until a desirable contour is found and a seed element selected.
  • the selection of a starting point to render an internal surface of interest can be done using computer selection means 34.
  • a threshold value supplied by the operator is used in sequentially scanning the data array for elements that cross the threshold value.
  • the first pixel having the same value as the threshold value corresponding to the surface of interest is selected as a starting point and its memory address is passed to processor 100 which generates an internal surface of interest.
  • processor 100 After completing the given surface, as indicated by absence of face-adjacent voxels that share a face intersect with voxels being intersected by the surface, processor 100 sends a signal back to selection means 34 prompting for selection of the next starting point.
  • selection means 34 Since all voxels visited by the processor 100 are "flagged" to indicate that they have been processed, selection means 34 only identifies the next unprocessed starting point on a surface of interest. Its address is passed to the processor for rendering the next surface structure. The process continues until no voxel within the stored data array having a vertex value corresponding to the desired threshold is left without being flagged as processed. This provides an indication that all disjoined surfaces that satisfy the threshold conditions have been processed.
  • Storage means 60 is used by the SpiderWeb processor to store data about the interconnection structure of the processed isosurface. In accordance with the processing method, this may include information about the visited voxel vertices, edges and faces, as well as the computed locations of edge hits and articulating points, and information about the interconnection of different elements into individually addressable data structures, such as the approximation triangles spun from each face intersect. Also stored in storage means 60 are the so-called shared normal vectors associated with the articulating points and voxel edge hit locations. The shared normal vectors are computed in processor 100 as the average of the set of unit normal vectors for each triangle having a vertex at the respective hit or articulating point.
  • the display of the generated approximation surfaces is an important element of the proposed invention. While both SpiderWeb A and B methods in general exhibit insensitivity to noisy data which is results in a better display, they both generate a large number of triangles, thereby increasing the computational burden of the system. In a number of cases, as discussed earlier, such triangles are in fact co-planar and in principle need not be processed separately. The specific data structure used in the methods may, however, be used to significantly reduce the computational burden at the display stage.
  • the T-mesh run illustrated in the figure is a contiguous strip of triangles, ordered in a "zig-zag” fashion.
  • a set of three initial triangle vertices, labeled 0, 1, 2 is first drawn.
  • the triangle 0 is drawn, using the connecting edges between vertices as delimiters. Since the edge (1-2) is shared by triangles 0 and l, when rendering the latter triangle, one need only specify one more point (vertex 3) .
  • vertex 3 is specified, triangle 1 can be drawn.
  • vertex 4 builds triangle 2, by only specifying one additional point. The process continues, until all triangles of the T-mesh run are drawn. At each step only one additional triangle vertex is specified and the display processor never needs to draw an edge twice.
  • the SpiderWeb method of surface rendering results in a data structure especially suitable for display using the T-run approach illustrated in Fig. 13. Since information about the shared triangle edges is already processed and stored, the appropriate T-mesh sequence can correspondingly be defined by traversing 5 voxel faces and considering rows of triangles as a group.
  • a Silicon Graphics VGX or VTX hardware can be used to draw surfaces at the rate of more than 1.6 x 10 6 T-mesh triangles per second.
  • the hit location is taken to a 5 standard position (a preferred choice is midway between the edge vertices) .
  • the independent angles between surface triangles are thus reduced to a relatively small number of cases, which simplifies the display at the expense of a somewhat "blocky" o appearance of the surface.
  • Figs. 15 and 16 illustrate the difference between a fully processed, and a blocky approximation surfaces.
  • the output structure of the SpiderWeb processor 100 can be described as 5 being in a vector form, as a sequence of triangle vertex points.
  • the data supplied to display processor 70 comprises a list of such vectors in three dimensions.
  • the list also includes the shared normal vector, 0 associated with each triangle vertex.
  • Such normal vectors are computed at steps 260 and 370 in Figs. 9 and 12, respectively.
  • a shaded, rasterized image is produced by the display processor 70 and supplied to 5 display 90, in which each triangular tile is shaded on a scale, dependent on the relation between the viewing angle selected from the control 50 and the normal vector to the approximation tile.
  • the unit normal of all triangles sharing a common vertex is averaged together at the vertex to produce a shared, or vertex normal.
  • Each triangle normal is computed as the vector cross product of two defining edges.
  • Figure 17 presents an illustration of the concepts.
  • the graphics processor 70 can then use the shared normals, to grade the calculated intensity over the triangle to create a smoothly shaded image having a natural appearance. Different coloring can also be used in the processor to provide better visual separation of the individual surface tiles, and of different objects or structures.
  • FIGS. 18 and 19 show examples of medical imaging surfaces, rendered in accordance with the method of the present invention.

Abstract

A system and method for surface rendering and for the display and analysis of arbitrary structures within the interior of a solid body is presented. The method uses as a starting point a voxel element (200) which is known to be crossed by a surface of interest (220). The surface is rendered in a topologically consistent way by connecting encountered surface hits (230) at each voxel face. A pair of voxel face hits, surrounding a vertex value exceeding the surface value are used to form a triangle, spun around a third articulating point in the interior of the voxel (240). Voxel faces sharing the surface hits are sequentially examined, and a set of approximation triangles is formed, to build a triangularization grid (250). The rendered surfaces are interactively displayed (280) and analyzed for global structure properties.

Description

SYSTEM AND METHOD FOR SURFACE RENDERING OF INTERNAL STRUCTURES WITHIN THE INTERIOR OF A SOLID OBJECT
FIELD OF THE INVENTION This invention generally relates to a system and method for surface rendering and more particularly to a system and method for the display and analysis of arbitrary three dimensional surface structures.
BACKGROUND OF THE INVENTION In many practical applications it is required to obtain and analyze three-dimensional (3D) arrays of data representing one or more physical properties at regular grid positions within the interior of solid bodies. Typically, the data is obtained by methods such as computer axial tomography (CAT) , x-ray scanning, nuclear magnetic resonance (N R) imaging systems, or by other non-intrusive mechanisms such as ultrasound scanning, positron emission tomography (PET) , emission computer tomography (ECT) and multimodality imaging (MMI) . Each of these techniques produces a planar, grid-like array of values for a continuous sequence of planes or slices of the solid object, thus providing a three-dimensional array of values. A typical application is in medicine where the solid object can be a human body or a part of it and the objective is to detect and model certain internal structures for diagnostic and surgical planning purposes. The methods to which the present invention relates are, however, equally applicable to other natural or artificial bodies and fields of study.
Examples of the physical properties measured at regular 3D positions include the coefficient of x-ray absorption in the CAT case or, for NMR imaging, the spin-spin or the spin-lattice relaxation parameters. In each case, the measured physical values reflect variations in the composition, density or structural characteristics of the underlying physical objects, providing knowledge about their internal structure. The 3D data array typically corresponds to a plurality of measurements of the given physical property, taken at regular intervals along an appropriately selected system of 3D (X, Y, Z) physical coordinates. Each measurement V^ of the physical property is thus associated with a specific array index position in a cubic lattice within the volume in study. A cubically adjacent set of eight such positions defines an unit cubic volume of the body under study and is known in the art as a voxel.
It is further known in the prior art to utilize such three-dimensional arrays of interior physical values to generate visual images of the structures within the body on a display system. In order to display two-dimensional (2D) images of such 3D interior structures, however, it is necessary to locate the position of the surface of such structure within the array of physical values. This is accomplished by comparing the array values to a threshold value or to a range of threshold values which correspond to the physical property values normally associated with that surface.
For example, in the case of the human body, bones, brain or any other tissues can be characterized by the known range of their density values which is compared to the measured and stored 3D array values in order to separate out a surface model of their structure. The displayed visual images can next be used to assist the medical practitioner in providing a diagnosis or in the planning of complex surgical procedures.
Once the location of the surface of an internal structure is determined, this surface is typically illuminated and shaded so as to give the human eye the impression of a natural shape and spatial disposition of the surface when it is displayed. Normally, to provide such shading the angular direction of a vector orthogonal to the surface at each point on the surface (the surface normal) is compared to the viewing angle of the observer. The intensity of shading can be adjusted to be proportional to the difference between these angles. There are several problems associated with the surface rendering of structures in the interior of a solid body. First, it is necessary to accurately determine which of the voxel elements of the 3D array have been crossed by the surface of interest and at which physical location. Next, it is important to provide a mechanism for displaying the internal structure of a three-dimensional body on the two- dimensional computer or TV screen. Another important problem is that data acquired from physical data collections systems is often noisy and is not capable of being directly displayed prior to some kind of a filtering operation.
Several solutions have been proposed in the past to address these problems. An approach that is frequently used is the one described in U.S. Pat. No.
4,710,876 to Cline et al. on the so-called "marching cubes". In this method, the surface segment intersecting a voxel is approximated by one of a limited number of standardized plane polygonal surfaces intersecting the voxel. One particular standardized surface tiling is selected from a library of available surfaces by a vector which represents the binary differences between the threshold value for the given surface of interest and the eight voxel vertex values. The surface-to-voxel intersection coordinates, as well as the vector normal to the surface for each such standardized polygonal surface tiling can be obtained by table look-up techniques. The final surface is assembled as a mosaic using all the standardized triangular tiles.
The marching cubes method has found significant acceptance in the industry because it is relatively simple to implement and is fast. Several modifications have been proposed since 1987. For instance, in U.S. Pat. No. 4,729,098 granted in March 1988 to Cline et al. a variation of the marching cubes methods is proposed, in which non-linear interpolation is used to more accurately locate the coordinates of the standardized polygons. Another known method for approximating the surface of an internal structure is the so-called "dividing cubes" method, disclosed in U.S. Pat. No. 4,719,585 to Cline et al. In this method, the values at the vertices of the voxel are used to interpolate in three dimensions values at regularly positioned intra-voxel sub-grid locations. The interpolated sub- grid values can then be used to locate the surface position and calculate the normal vector more accurately. In order to distinguish between different internal structures having the same or similar physical property and thus having close or similar values on the three-dimensional grid, U.S. Pat. No. 4,751,543 to Lorensen discloses a technique for labeling surfaces with similar property values and using adjacency criteria in respect to a "seed location" in the particular structure of interest to segregate the desired surface from all of the labeled surfaces.
U.S. Pat. No. 4,791,567 to Cline et al. discloses another technique of segregating similar internal structures by determining connectivity from adjacency information. U.S. Pat. No. 4,821,213 to Cline et al. discloses yet another technique for differentiating internal structures in which a linear path is made through the data array to locate and label all different structures along the scan line by counting structure interfaces.
In the methods described above, voxels within the solid object under study are typically traversed sequentially and progressively, so that every voxel within the data space is checked for a crossing by the surface of interest, and is subsequently assigned to belong either to the outside or to the inside of the surface. Because of the large volumes of data to be processed, the methods are computationally intensive and relatively slow. An improvement in this respect is described in U.S. Pat. No. 5,166,876 to Cline, issued in November 1992. The disclosed method significantly speeds up the rendering of a surface within a boundary by only processing voxels on the surface. This approach saves significant computation time but still uses the same case table and the limited set of approximating surfaces which were described in the original marching cubes patent. As noted above, the marching cubes surface rendering method is widely accepted but has certain deficiencies. These deficiencies are generally related to the fact that the fifteen canonical cases used in the surface segments' approximation do not cover the entire spectrum of possible tiling combinations. This results in occasional approximation errors, occurring at lower levels of resolution, specifically if a voxel face is intersected by the surface twice resulting in four crossings ("hits") with the voxel edge. The errors manifest themselves in incorrect surface rendering, visually appearing as unexpected spots, tears and missing tile elements. The problems are aggravated if the collection data is noisy and/or with low contrast. A number of corrections to the marching cube algorithm have been offered in the past. Notable is a paper by Neilson and Hamann, "The Asymptotic Decider: Resolving the Ambiguity in Marching Cubes," in the Proceedings of Visualization 1991, pp. 83-91. This and other publications show that the case table approach produces inconsistent results in some cases. Thus, no topologically correct algorithm for isosurface tiling has been offered so far. Accordingly, reliable measurements of internal structure parameters, such as surface area and volume, are also not available.
The display of 3D graphic images has been principally driven by the goals and directions of computer aided design (CAD) and computer aided manufacturing (CAM) systems. Systems accepting data in several basic formats have been developed for displaying solid bodies and manipulating images in various ways to create models of parts being manufactured, or internal organs being monitored. The present systems for surface rendering and display are, however, still incapable of meeting the requirements for a fast and truly interactive display of continuously monitored 3D data. The complexity of the task can be better appreciated by noting that a typical surface rendering operation involves processing on the order of 512 x 512 x 512 pixels. In the presence of noise the task becomes even more complicated because of tbe necessity for additional noise filtering, smoothing and data interpolation, which takes processing time and often results in the loss of important structural details.
It is thus desirable to provide a method for rendering topologically correct surface approximations of arbitrary internal structures within a volume of interest and a system for the display and analysis of such surface structures. It is also desirable to provide high quality, interactive display of the surface structures of interest from a set of collected data and minimize the effects of possible noise interference. Additionally, it is desirable to provide an efficient implementation of the method which facilitates the use of specialized hardware for real time data processing.
SUMMARY OF THE INVENTION Accordingly, it is an object of the present invention to provide a system and method for topologically correct surface rendering of internal structures within a volume of interest.
It is another object of the present invention to provide a system and method for the display of three- dimensional information. It is also an object of the present invention to provide a system and method which can be implemented in specialized electronic hardware.
It is another object of the present invention to provide a graphics imaging system capable of being interactively used, and providing high quality images with textural shading and other visual clues to the user.
It is yet another object of the present invention to provide medical practitioners with the abilities to emulate surgical procedures graphically prior to undertaking invasive patient treatments.
Additionally, it is an object of the present invention to provide a plurality of three-dimensional surface views from a single set of collected data.
A further object of the present invention is to provide a system for analysis of internal structures within a 3D solid body.
Lastly, but not exclusively it is an object of the present invention to minimize the effects of possible noise interference.
In accordance with a preferred embodiment of the present invention, a novel method for rendering three- dimensional surface structures comprises the following steps. Signal patterns representing the value of at least one physical property associated with regularly spaced grid locations in a three dimensional body are stored in a memory means. The regularly spaced grid locations define voxel elements within the body, such that each voxel is associated with eight vertex values and six voxel faces. The signal patterns that represent values lying within a specified range associated with the surface structure of interest are next selected for further examination. An initial grid location which corresponds to a position within or on the surface structure is first selected and the eight vertex values associated with a voxel at the initial grid location are selectively retrieved. A threshold value is supplied to determine which voxel vertex values are inside and which voxel vertex values are outside of the surface structure. For each voxel face having at least one voxel vertex inside and at least one voxel vertex outside of the surface structure, there must be at least one pair of intersects between the edges of the voxel face and the surface structure. The coordinate values of these intersects between the voxel face and the surface structure are then computed. Preferably, the computation of these intersect coordinate values is done by a linear interpolation of the voxel vertex values. Other methods of computation can similarly be used.
A triangular approximation to the surface structure is then formed using the pair of intersects between the edges of the voxel face and the surface structure as two vertices of the triangle with a third vertex in the interior of the voxel. The step of forming triangles is repeated for each pair of voxel face intersects on the six faces of the voxel. The physical location of the third vertex of each triangle is identical and is computed as the centroid of all interpolated intersect coordinate values for the voxel.
The method further comprises the steps of identifying all adjacent voxels sharing a face intersect with the surface structure. Next, for each identified adjacent voxel, the steps of retrieving the eight vertex values, determining the inside and outside vertex values, computing the coordinates of the voxel face intersects and forming the triangles is repeatedly applied until all face-adjacent voxels within the three-dimensional body having a face intersect with the surface structure are processed. In the step of computing the intersect coordinates, only values that have not previously been computed are processed. The proposed method of surface rendering further comprises the step of displaying a two-dimensional image of the approximated surface structure. Preferably, the surface structure is displayed by utilizing only the voxels identified by the processor as being crossed by the surface structure.
In a second preferred embodiment of the present invention, a method for generating of an internal surface structure within a three-dimensional body comprises the following steps. Three-dimensional signal patterns representing the value of at least one physical property associated with regularly spaced grid locations in a three dimensional body are stored in memory means. The regularly spaced grid locations define voxel elements within the body, such that each voxel is associated with vertex locations and voxel faces, each two intersecting voxel faces defining a voxel edge. Each voxel edge passes through two voxel vertices.
The signal patterns which represent values lying within a specified range associated with the surface structure of interest are next selected. Two voxel vertices of a voxel edge for one of said stored voxels are then retrieved; and it is determined which voxel vertex of the selected voxel edge is inside and which voxel vertex of the selected edge is outside of the surface structure.
For each selected voxel edge having at least one voxel vertex inside and at least one voxel vertex outside of the surface structure, a coordinate value of a first hit location where the surface structure crosses the voxel edge is next computed. Preferably, the computation of the intersect coordinate values is done by a linear interpolation of the voxel vertex values but other computation methods can be used as well.
Next, each of the two intersecting voxel faces that define the voxel edge where a hit was identified is examined for a second hit location which defines an intersect between the surface structure and the voxel face.
For each intersect of a voxel face with the surface of interest, a triangle with a vertex in the interior of the voxel is formed to obtain an approximation to the surface structure. The location of the third vertex of each approximating triangle is computed as the centroid of the voxel face intersects such that the voxel face intersects form a closed contour around the inside voxel vertex.
The method further comprises the steps of identifying all adjacent voxels that share a face intersect with the surface structure. Next, for each identified adjacent voxel, the steps of selectively retrieving two voxel vertices of a voxel edge, determining the inside and outside vertex values, computing coordinates of the first hit locations, examining the voxel faces defining the voxel edges, and forming the approximation triangles are repeatedly applied until all face-adjacent voxels within the three-dimensional body having a face intersect with the surface structure are processed. In the step of computing the intersect coordinates, only values that have not previously been computed are processed.
In accordance with the present invention, the method further comprises the step of displaying a two- dimensional image of the approximated surface structure. Preferably, the surface structure is displayed by utilizing only the voxels identified by the processor as being crossed by the surface structure.
In accordance with a preferred embodiment of the present invention, the proposed method for 3D surface rendering is implemented on a system for analysis and display of three-dimensional surface structures which comprises the following elements. First, storage means are provided for storing 3D signal patterns representing the value of at least one physical property associated with a 3D body at regularly spaced grid locations within the body. The signal patterns are provided by a three-dimensional raw data source, such as MRI, CAT, x-ray, ultrasound or others. A set of 3D-adjacent grid locations define a volume element referred to as voxel having associated with it voxel edges and voxel faces.
The system also comprises means for defining a plane of pixel positions having preselected values and corresponding to one plane of the spaced grid locations. Also provided is means for determining in the defined pixel plane, a pixel associated with a voxel that is intersected by a surface structure of interest within said body. The system further comprises processing means for computing coordinate values of the intersection between the surface structure of interest and the intersected voxel edges and for computing the coordinate values of at least one articulating location positioned in the interior of the intersected voxel. Second storage means for storing the computed coordinate values along with connectivity information about joining the computed coordinate values into a set of triangles is also provided. These triangles approximate the intersection between the surface structure of interest and the voxel. Each triangle has one edge on a voxel face and one vertex at an articulating location in the interior of the intersected voxel. Preferably, the processor means, further comprises generating means for generating normal vectors associated with the computed coordinate values. The computed normal vectors are similarly stored in the second storage means.
In a preferred embodiment of the present invention, the processing means of the system further comprises means for identifying which of the voxels sharing a voxel face with the intersected voxel are intersected by the surface structure of interest. To provide for a better data management of the system, the processing means further comprises means for generating, for each voxel that is intersected, a plurality of index values representing information about the shared voxel face intersects. These index values are also stored in the second storage means and are used to indicate to the processor which voxels and voxel faces have been processed.
The processing means of the system comprises means for simultaneously comparing all vertex values of an intersected voxel to a threshold value. This comparison is done independently for each voxel face thus enabling a parallelizable, high speed implementation of the proposed system.
In order to facilitate the display of the 3D surface structure of interest, the system of the present invention further comprises display processor means which receives data from the second storage means and converts it into a display format. The display format data is supplied to means for interactively displaying surface structures of interest. In accordance with another embodiment of the present invention, the planar arrangement of the regularly spaced grid-locations is other than a rectangular array. For example, the planar arrangement of the grid locations may have a hexagonal, or other pattern suitable for processing shape.
The system and method for surface rendering and display of the present invention have a number of practical applications. For example, the internal data structure representation allows the user to simultaneously process and display two or more isosurfaces. As a result, the system provides a mechanism to count 3D disjoint structures within the interior of a body under study by enumerating the independent structures being processed.
The system and method of the present invention can in addition be used to compute and display cross- sections of disjoint surface structures. Such surface structures need not only be internal structures of the object under study. While the description of the proposed method is presented in terms of the measured values of a given physical property of the object, other types of surfaces can similarly be computed and displayed. For instance, this can include isogradient and isocurvature surfaces giving information about the locations within the body having identical rate of change in the density values or curvature values, respectively. The cross-section of such surfaces with the structure surface can be used to automatically compute cardinal landmarks on the object surface.
In a different application, the system of the present invention can be used to compute global structure parameters such as the total surface area, and the structure volume. Information necessary to compute these quantities, including the coordinates of a set of physical points which belong to the surface is stored and available for each processed structure. The resulting parameter estimates can be used to define characteristics of the object under study, such as volume occupancy, average structure size and others, which characteristics have a wide range of applications, especially in 3D pattern classification. Lastly, but not exclusively, the system of the present invention can be used to observe time changes in the interior structures of a given solid object by knitting together "snap-shots" of the surface of the structure and displaying these snap-shots in a sequence. An exemplary application of such system use is the monitoring of a patient's heart beat pattern.
BRIEF DESCRIPTION OF THE DRAWINGS These and other advantages, features and objects of the present invention are described in the following description of the preferred embodiment in which:
Fig. 1 is a detail of surface tiling errors in a prior art surface rendering method.
Fig. 2 is a block-diagram of a system for carrying out the method of the present invention. Fig. 3 is an isocontour plot of an MRI axial slice of a brain used for plane selection preview. Fig. 4 is a three dimensional voxel with elements defined.
Figs. 5(A) and (B) are perspective views of a voxel neighborhood and voxel adjacencies for three dimensional voxel elements, respectively. Fig. 6 shows intersection of a voxel by surface structures. Fig. 7 illustrates voxel face hits assumptions. Fig. 8 shows ambiguity in voxel face contour lines. Fig. 9 is a flow-chart of the SpiderWeb A method. Fig. 10 illustrates spinning approximation triangles using an articulation point.
Fig. 11 illustrates spinning approximation triangles for a voxel face contour ambiguity using SpiderWeb A.
Fig. 12 is a flow-chart of the SpiderWeb B method.
Fig. 13 shows a two-dimensional adaptive sampling grid. Fig. 14 shows triangle vertices drawn in a T-mesh order.
Fig. 15 is an approximating surface computed with full precision.
Fig. 16 shows the approximating surface of Fig. 16 using midway hit location speed processing.
Fig. 17 shows triangle face and vertex normal vectors.
Figs. 18 and 19 illustrate the medical imaging surfaces rendered in accordance with the method of the present invention.
DETAILED DESCRIPTION OF THE INVENTION Referring now to Figure 1 a surface detail is shown which illustrates tiling errors in a prior art surface rendering method. The errors are observed as bright or very dark spots on the surface, indicating a tear in the represented surface or a missing approximation tile. The present invention is designed to remove such errors by presenting a topologically consistent method of building a triangularization grid which represents an arbitrary three-dimensional structure.
Illustratively, the invention is practiced on the raw data obtained from non-intrusive measurements of at least one physical property within a planar slice through a body under study. A typical source of the raw data is an MRI or CT system, an x-ray, ultrasound or another type of a scanner, capable of converting measurements of some physical property of the body, at regular grid-like locations within the planar slice, into an array of discrete values. Other non-intrusive data collection mechanisms such as positron emission tomography (PET) , emission computer tomography (ECT) and multimodality imaging (MMI) systems can also be used as a data source.
The operation of the data source may be carried out for a sequence of slices of the solid object, thus providing a three-dimensional array of values representative of the spatial structure of the object. Alternatively, it is possible to collect data over various periods of time, for the purpose of measuring temporal changes in the body under study.
The raw data is stored in a first storage means, preferably indexed by a slice number. The storage means preferably has a fast access time and may comprise magnetic storage, random access memory or other known in the art storage media. Advantageously, the storage means also comprises means for converting the raw data into an application specific internal data structure, in order to take advantage of the specific hardware configuration used. For example, the internal data structure provides the capability of individually addressing and retrieving individual pixel elements as well as sets of spatially adjacent pixel configurations, referred to as voxels.
In accordance with a preferred embodiment of the present invention, the first storage means also provides the capability of indexing each vertex, each voxel edge and voxel face, in order to facilitate the computation process.
The indexing of the voxels will be better understood from Figs. 4-6. Fig. 4 shows a voxel element consisting of a cubic arrangement of eight data pixels and presents element definitions consistently used in the following description. Each pixel element is associated with a triplet of indexes (i,j,k) used to determine the location of the pixel in a three-dimensional indexing system. By convention, the lower left vertex of the front face is selected as a reference point. The indexes of the adjacent pixels to the left, to the top and towards the back plane are shown to increase by the length s, which corresponds to the number of display pixels inscribed within two consecutive data points. It is to be understood .that Fig. 4 only serves as an illustration and other possible conventions are possible, as long as they consistently applied. Thus, a different vertex of a voxel element may be selected as a reference point, and the direction of index increase may be selected opposite to the one shown in Fig. 4. Furthermore, the length s may correspond to substantially different physical distances depending on the particular application. In addition, even though in Fig. 4 the length is identically labeled in all three index directions, in practice this may not necessarily hold true. The reason is that the data sampling step in each coordinate direction depends on the corresponding spatial resolution of the raw data source in the specific index direction. In accordance with the present invention, voxel elements, such as the vertices, edges and voxel faces, illustrated in Fig. 4, can be referenced by a data processor in groups or individually. Label information is associated with each voxel element, so that vertex values connected along edges of the voxel, or voxel faces are marked as being connected. Similarly, other data structures which are used in a group during the processing are marked as being connected. A second storage means associated with the data processor comprises a memory which is allocated for each data structure, so that the processor is capable of determining whether or not a particular element of the data structure has been previously processed. For instance, each time a specific voxel edge is retrieved, a mark is left indicating for all associated data structures that this element has been processed. Thus, when the next element of a connected data structure is visited, the processor need not revisit the original voxel edge. Fig. 5 (A) is a perspective view of a voxel neighborhood. Fig. 5 (B) depicts possible combinations for voxel pairs, which define a vertex adjacency, an edge adjacency and a face adjacency, depending on whether only a vertex, an edge or a voxel face, respectively, are shared. In reference to Fig. 6, a voxel element is shown intersected by two separate isosurfaces. Voxel vertices, having values above an isosurface threshold are termed "inside" vertices, while voxel vertices having values less than the isosurface threshold are termed "outside" vertices. The intersection of an isosurface of interest with a voxel edge is referred to as a "hit". Six such hits, indicated as H0-H5 are shown in Fig. 6. The intersection between the surface of interest and a voxel face is referred to as a "face intersect" and, as seen in Fig. 6, connects two voxel edge hits.
As previously noted, the first storage means comprises means for converting the raw data into a specific internal data structure, which allows all voxel elements of the data array to be consistently indexed by the data processor for identification purposes. For example, the face intersect H0-H1 on the left face of the voxel in Fig. 6 is shared with the adjacent left voxel and, in fact, occupies the same physical location. By using a consistent indexing scheme, as well known in the computer graphics art, the data processor need only compute the hit locations once, and use their values when processing all corresponding edge adjacent voxels, as defined in Fig. 5(b).
In accordance with the present invention, the data processor implements a topologically consistent method of connecting the isosurface hits on voxel edges into a continuous triangularization grid representing the surface of interest for analysis and display. The method requires a starting point, or a "seed" pixel element on a data plane, selected by a plane selecting means. Using only local information about the adjacent pixels, the method "spins" a web of triangles, closely approximating the surface of interest. The processing sequence and visual appearance of the rendered surface resemble the construction of spider web, giving the process its "Spiderweb" name.
Some assumptions about the topology of the approximated surfaces used in the surface rendering method of the present invention are illustrated in reference to Figures 7 and 8. In Fig. 7 depicts a voxel face with vertices V0-V3, edges E0-E3 and a clockwise scanning direction, as indicated by the arrow. A hit HI is shown on voxel edge El. In accordance with the present invention, it is first assumed that the crossing isosurface of interest is sufficiently smooth. This assumption imposes a limitation on the surface curvature.
Secondly, an isosurface is assumed to have no singularities, breaks, or holes, so that given a hit on one edge, there will be at least one more hit on the same voxel face. As a result of the these assumptions, in accordance with the present invention hits on a voxel face must come in pairs. Referring back to Fig. 7, this implies that there must be at least one additional hit on the voxel face. If no i such hit, or an odd number of hits is found, the processor will indicate an error.
Figure 8 illustrates another possible situation, when all four edges of a voxel face have isosurface hits. As shown, this is only the case when diagonally opposite voxel vertices are both "inside," or both "outside" the surface of interest, and the two diagonal pairs are not on the same side of the surface. Lacking global information about the isosurfaces, the exact direction of surface intersects cannot be determined. This ambiguity results in four possible cases of face intersects, indicated as cases A, B, C and D in Fig. 8. Cases B and D will be re¬ considered in connection with the SpiderWeb A, and SpiderWeb B processing methods described next.
SpiderWeb A
Figure 9 is a flow-chart of a first embodiment of the surface rendering method in accordance with the present invention, referred to as SpiderWeb A. In this method, no decision about the exact direction of the face intersects in the case of four hits is made. This corresponds to case D in Fig. 8, where the four hits are connected, thus minimizing the probability that an error is made. The algorithmic steps in this case are as follows.
In step 200 of the method, the pixel elements corresponding to one voxel are simultaneously retrieved from the data storage. The vertex elements f the voxel are next compared in step 210 to a preselected threshold which corresponds to the surface of interest, and each vertex is assigned as either inside or outside with respect to the surface.
If an intersection of the surface with the voxel is not detected, the method proceeds to retrieve a face adjacent voxel in accordance with predefined voxel sequencing rules. If a surface intersection is detected at step 220, the locations of all edge hits of the voxel, which have not been previously processed, are computed. Preferably, the position of a hit along an edge is computed by using a linear interpolation of the values at the two vertices defining the edge. For example, if a voxel edge with a hit has an outside vertex at position Vx0, with a value of V0, and an inside vertex at position Vxl with a value of V,, and the threshold value is T, SO that V0 < T < V,, the location of the edge hit Ex can be calculated using the formula:
Ex = V^ + (Vxl - Vrf) (T - V0)/(V, - V0)
Naturally, other methods of estimating the position of the hit, such as non-linear interpolation, may be used, as is well known in the art. The next step 240 involves the calculation of an Articulating Point (AP) which is located entirely in the interior of the voxel. The articulating point in the SpiderWeb A method is computed as the mean (centroid) of the system of all connected hits in a voxel.
The articulating point is used in the next step 250 as the third vertex of a set of surface approximation triangles, each triangle being defined by the articulating point and each of the possible different pairs of the connected hits in the voxel. As can be appreciated by those skilled in the art, if only two edge hits are found for each voxel face, as for example shown in Fig. 6, the articulating point will lie on the plane defined by the hits; and the approximation triangles will be co-planar. An illustration of the SpiderWeb A two-face-hits approximation surface is shown in Fig. 10 where the articulating point is located near the center of the voxel at the centroid of the four hits and four approximation triangles are defined, each with the articulating point as one vertex and with two connected hits on one of the four faces as the other two vertices.
If, on the other hand, four hits are encountered in a voxel face, the resulting approximation surface triangles will form a pyramid, touching at its tip (the AP) another approximation surface. An example is shown in Fig. 11.
In step 260 of the method, triangle and vertex normals for each of the appr.oximation triangles are computed. The purpose of these normals will be discussed in more detail below in connection with the description of display processor 70 of the system of Fig. 2. In the following step 270, for each intersected voxel face, the face adjacent voxel is retrieved and processed. In accordance with the present invention, the indices of all voxel faces that have already been processed are marked so that processed voxel faces need not be considered again.
Finally, data from the SpiderWeb processor is forwarded to storage and display at step 270.
SpiderWeb A processing method has the disadvantage of creating "tunnel" voxels and "bubbles" in the case of a voxel face ambiguity. However, in the case of a noisy data, this property may actually reduce the visual impact of the noise, since its effect is localized to a particular voxel element and is not visible to the user.
SpiderWeb B
Figure 12 is a flow-chart of the second embodiment of the surface rendering method in accordance with the present invention, referred to as SpiderWeb B. In this method, a consistent decision is made about the exact direction of the face intersects in the case of four hits on a voxel face. The decision corresponds to case B in Fig. 8, where the four edge hits are connected in a way so as to surround the inside voxel vertex. The algorithmic steps in this case are as follows. In step 300 of the method, a starting or seed vertex is selected. In the next step 310, the three edges associated with the seed vertex are scanned so as to retrieve the vertex elements at the opposite ends of those edges, and in the following step 320 those vertex elements are compared to the pre-selected threshold. The hit locations are calculated by means of a linear interpolation, as discussed above, or estimated in any other suitable manner, as known in the art.
Triangles are formed in SpiderWeb B in the following four steps which are performed for each edge that is hit. First, in step 340, the two voxel faces that define the edge that has been hit are scanned for additional hits. This is done beginning with the inside vertex of the edge that has been hit and proceeding from vertex to vertex around the edges of the first voxel face and then around the edges of the second voxel face. Hit locations are calculated by linear interpolation as discussed above. Since face intersects on a voxel face occur in pairs, the scanning of each face stops when the first additional hit is located. The processor then forms a triangle having one vertex at the intersection point on the edge that was originally hit, a second vertex at the location of the first additional hit and a third vertex at an articulating point with an initially undetermined physical location inside the voxel. If the foregoing procedure identifies any single hits on another face of the voxel, that face is likewise scanned to locate the second hit; and the two hits on that face and an articulating point at an undefined location within the voxel constitute another surface approximation triangle. If a pair of hits are identified on another face of the voxel, these hits and an articulating point likewise constitute another surface approximation triangle. This process .^continues until all surfaces of the voxel are 5 accounted for. The exact location of the articulating point for all triangles of the voxel is then computed as the centroid of all edge hits, and information about the triangles is stored. Triangle and vertex normals, to be discussed next, are computed at step 0 370.
For each voxel edge with a surface hit, the process repeats itself in the face adjacent voxels. In particular, in step 380 a set of edges from the face adjacent voxel is examined until all possible 5 intersections of the isosurface of interest are processed. Finally, in step 390 the processor exits to storage means 60 and the display processor 70.
Unlike the SpiderWeb A method, this embodiment of the present invention does not necessarily involve o simultaneous processing of entire voxels. Instead, after the initial seed pixel is specified, the proposed method tracks along the isosurface of interest in accordance with the adopted indexing scheme, only visiting voxel edges and faces with hits. 5 Both SpiderWeb A and B processing methods reflect the uncertainty of the voxel face contour line ambiguity, illustrated in Fig. 8. In the embodiment of the present invention, using SpiderWeb A, no assumption is made about the exact direction of the 0 face intersects, resulting in a three dimensional "bubble". In the second embodiment of the present invention, using SpiderWeb B, a systematic choice to "surround inside vertex" is made based on physical considerations. Yet another possibility to avoid the 5 face intersect ambiguity is to increase the resolution of the space sampling grid until no voxel face has four hits per face. Thus, in accordance with another embodiment of the present invention the resolution of the sampling grid is adaptively refined. Figure 13 shows a two-dimensional adaptive sampling in which the upper left corner is illustratively subdivided until adequate resolution is achieved.
The rectangular array of sampling elements considered thus far is only illustrative of the different sampling grids which can be used in accordance with the present invention. Well known in the literature is for instance the two-dimensional hexagonal sampling grid that can reduce the number of required samples in a plane by about 13.4% (see, for example Dudgeon et al., "Multidimensional Digital
Signal Processing," 1984, pg. 39-47). Other sampling grids can also be used, as known in the art.
It is important to note that while the processor computes a large number of triangles, information for which needs to be stored, the actual computation process can be significantly sped up due to the fact that separate voxel faces can be processed independently. Thus, the proposed method can be implemented in a parallel architecture resulting in a significant improvement in terms of processing speed. As shown in Fig. 2, illustrative apparatus for the practice of the invention comprises a first storage means 20, a plane selection means 30, control unit 50, a Spiderweb processor 100, a second storage means 60, a display processor 70 and a structure analysis means 80.
Plane selection means 30 operates under the control of an operator from control unit 50 to select data from individual data planes stored in data storage 20. The operator may preview data planes corresponding to a particular physical location of interest in the body under study. Since structures having different physical properties exhibit different ranges of measured values, the operator can isolate a specific structure of interest by selecting a particular threshold or a range of threshold values. For example, the operator could isolate bones from the surrounding soft tissue in medical imaging.
Figure 3 illustrates an isocontour plot of an MRI axial slice of the brain used for a plane selection preview in accordance with a preferred embodiment of the present invention. In this embodiment, using selection means 32 in Fig. 2, the operator can choose a specific plane of interest and display it on a monitor. A preselect threshold value is supplied to display processor 70 and is used to display the locations (the contour lines) at which the selected data plane is crossed by a surface structure having the same physical property value. The operator can then identify a particular contour line and specify the coordinates of a desired starting point ("seed") . The coordinates of the seed are next supplied to processor 100 to initiate the rendering of the desired surface structure. If no contour line of the surface of interest is identified, the operator can interactively specify a different data plane or modify the threshold value until a desirable contour is found and a seed element selected.
In an alternative embodiment of the present invention, the selection of a starting point to render an internal surface of interest can be done using computer selection means 34. In this embodiment, a threshold value supplied by the operator is used in sequentially scanning the data array for elements that cross the threshold value. The first pixel having the same value as the threshold value corresponding to the surface of interest is selected as a starting point and its memory address is passed to processor 100 which generates an internal surface of interest. After completing the given surface, as indicated by absence of face-adjacent voxels that share a face intersect with voxels being intersected by the surface, processor 100 sends a signal back to selection means 34 prompting for selection of the next starting point. Since all voxels visited by the processor 100 are "flagged" to indicate that they have been processed, selection means 34 only identifies the next unprocessed starting point on a surface of interest. Its address is passed to the processor for rendering the next surface structure. The process continues until no voxel within the stored data array having a vertex value corresponding to the desired threshold is left without being flagged as processed. This provides an indication that all disjoined surfaces that satisfy the threshold conditions have been processed.
Storage means 60 is used by the SpiderWeb processor to store data about the interconnection structure of the processed isosurface. In accordance with the processing method, this may include information about the visited voxel vertices, edges and faces, as well as the computed locations of edge hits and articulating points, and information about the interconnection of different elements into individually addressable data structures, such as the approximation triangles spun from each face intersect. Also stored in storage means 60 are the so-called shared normal vectors associated with the articulating points and voxel edge hit locations. The shared normal vectors are computed in processor 100 as the average of the set of unit normal vectors for each triangle having a vertex at the respective hit or articulating point.
The display of the generated approximation surfaces is an important element of the proposed invention. While both SpiderWeb A and B methods in general exhibit insensitivity to noisy data which is results in a better display, they both generate a large number of triangles, thereby increasing the computational burden of the system. In a number of cases, as discussed earlier, such triangles are in fact co-planar and in principle need not be processed separately. The specific data structure used in the methods may, however, be used to significantly reduce the computational burden at the display stage.
Reference is now made to Fig. 13 in which a "T- esh" ordering of approximation surface triangles is shown. The T-mesh run illustrated in the figure is a contiguous strip of triangles, ordered in a "zig-zag" fashion. A set of three initial triangle vertices, labeled 0, 1, 2 is first drawn. After the last triangle vertex 2 is specified, the triangle 0 is drawn, using the connecting edges between vertices as delimiters. Since the edge (1-2) is shared by triangles 0 and l, when rendering the latter triangle, one need only specify one more point (vertex 3) . Thus, after vertex 3 is specified, triangle 1 can be drawn. In the next step, vertex 4 builds triangle 2, by only specifying one additional point. The process continues, until all triangles of the T-mesh run are drawn. At each step only one additional triangle vertex is specified and the display processor never needs to draw an edge twice.
The SpiderWeb method of surface rendering results in a data structure especially suitable for display using the T-run approach illustrated in Fig. 13. Since information about the shared triangle edges is already processed and stored, the appropriate T-mesh sequence can correspondingly be defined by traversing 5 voxel faces and considering rows of triangles as a group. A Silicon Graphics VGX or VTX hardware can be used to draw surfaces at the rate of more than 1.6 x 106 T-mesh triangles per second.
In order to further reduce the computational 0 requirements of the SpiderWeb method it is possible to simplify the computation of the edge hit locations. For example, if an edge has a hit, instead of calculating its location exactly, in one embodiment of the present invention the hit location is taken to a 5 standard position (a preferred choice is midway between the edge vertices) . The independent angles between surface triangles are thus reduced to a relatively small number of cases, which simplifies the display at the expense of a somewhat "blocky" o appearance of the surface. Figs. 15 and 16 illustrate the difference between a fully processed, and a blocky approximation surfaces.
Referring again to Fig. 2, the output structure of the SpiderWeb processor 100 can be described as 5 being in a vector form, as a sequence of triangle vertex points. In one embodiment of the present invention the data supplied to display processor 70 comprises a list of such vectors in three dimensions. The list also includes the shared normal vector, 0 associated with each triangle vertex. Such normal vectors are computed at steps 260 and 370 in Figs. 9 and 12, respectively. In a preferred embodiment of the present invention, a shaded, rasterized image is produced by the display processor 70 and supplied to 5 display 90, in which each triangular tile is shaded on a scale, dependent on the relation between the viewing angle selected from the control 50 and the normal vector to the approximation tile.
Specifically, in order to produce smoothly shaded images, the unit normal of all triangles sharing a common vertex is averaged together at the vertex to produce a shared, or vertex normal. Each triangle normal is computed as the vector cross product of two defining edges. Figure 17 presents an illustration of the concepts. The graphics processor 70 can then use the shared normals, to grade the calculated intensity over the triangle to create a smoothly shaded image having a natural appearance. Different coloring can also be used in the processor to provide better visual separation of the individual surface tiles, and of different objects or structures.
Since the computation of the vertex normals in accordance with the present invention uses localized planar triangle information, as opposed to other known methods of computing surface normals which typically use the volume gradient case, when closely positioned isosurfaces need to be discriminated, the method of the present invention provides a better surface separation. Figures 18 and 19 show examples of medical imaging surfaces, rendered in accordance with the method of the present invention.
It is to be understood that the claimed invention is not to be limited to the description of the preferred embodiment but encompasses other modifications and alterations within the spirit of the inventive concept which scope is defined by the following claims.

Claims

What is claimed is:
1. A method for rendering three-dimensional surface structures, comprising the steps of:
a) storing three-dimensional signal patterns representing the value of at least one physical property associated with a three dimensional body at regularly spaced grid locations defining voxel elements within said body, each voxel having eight vertex values and six faces;
b) selectively retrieving the eight voxel vertex values of said physical property for one of said stored voxels;
c) determining which voxel vertex values are inside and which voxel vertex values are outside of the surface structure;
d) computing, for each voxel face having at least one voxel vertex inside and at least one voxel vertex outside of the surface structure, coordinate values of the intersect between the voxel face and the surface structure; and
e) forming, for each intersect of a voxel face with the surface structure, a triangle having a third vertex in the interior of the voxel to obtain an approximation of the surface structure.
2. The method of claim 1 wherein step c) further comprises the step of comparing each voxel vertex value with a predetermined threshold, the value of said threshold corresponding to the surface structure.
3. The method of claim 2 wherein at step d) the computation of the intersect coordinate values is done by a linear interpolation of the voxel vertex values.
4. The method of claim 3 wherein at step e) the physical location of the third vertex of each approximating triangle is identical and is computed as the centroid of all interpolated intersect coordinate values.
5. The method of claim 4 further comprising the steps of: f) identifying each adjacent voxel, sharing a face intersect with the surface structure; and g) recursively applying steps b) , c) , d) , e) and f) to each identified voxel until all face-adjacent voxels within the three-dimensional body having a face intersect with the surface structure are processed.
6. The method of claim 5 wherein in the step of recursively computing the intersect coordinates, only coordinate values that have not previously been computed are processed.
7. The method of claim 5 further comprising the step of displaying a two-dimensional image of the approximated surface structure.
8. The method of claim 7 wherein the surface structure is displayed by utilizing only the identified voxels.
9. A method for rendering of an internal surface structure within a three-dimensional body, comprising the steps of:
a) storing three-dimensional signal patterns representing the value of at least one physical property associated with a three dimensional body at regularly spaced grid locations defining voxel elements within said body, each voxel having associated with it vertex locations and voxel faces, each two intersecting voxel faces defining a voxel edge, said voxel edge passing through two voxel vertices;
b) selectively accessing two voxel vertices of a voxel edge for one of said stored voxels;
c) determining which voxel vertex of the selected voxel edge is inside and which voxel vertex of the selected edge is outside of the surface structure;
d) computing, for each selected voxel edge having at least one voxel vertex inside and at least one voxel vertex outside of the surface structure, a coordinate value of a first hit location where the surface structure crosses the voxel edge;
e) examining each of the two intersecting voxel faces defining the selected voxel edge for a second hit location defining an intersect between the surface structure and the voxel face, said intersect passing through said face, said intersect passing through said first and said second hit locations;
f) forming, for each intersect of a voxel face with the surface of interest, a triangle with a vertex in the interior of the voxel to obtain an approximation to the surface structure. 0
10. The method of claim 9 wherein step c) further comprises the step of comparing each voxel vertex value with a predetermined threshold, the value of said threshold corresponding to the surface structure. 5
11. The method of claim 10 wherein at step d) the computation of the intersect coordinate values is done by a linear interpolation of the voxel vertex values.
o 12. The method of claim 11 wherein at step e) the location of the third vertex of each approximating triangle is computed as the centroid of the interpolated intersect coordinate values forming a closed contour around the inside voxel vertex. 5
13. The method of claim 12 further comprising the steps of: g) identifying each adjacent voxel, sharing a face intersect with the surface structure; and 0 h) recursively applying steps b) , c) , d) , e) and f) to each identified voxel until all face-adjacent voxels within the three-dimensional body having a face intersect with the surface structure are processed.
5 14. The method of claim 13 wherein in the step of recursively computing the intersect coordinates, only coordinate values that have not previously been computed are processed.
15. The method of claim 13 further comprising the step of displaying a two-dimensional image of the approximated surface structure.
16. The method of claim 15 wherein the surface 0 structure is displayed by utilizing only the identified voxels.
17. The method of claim 9 wherein the planar arrangement of the regularly spaced grid-locations is 5 other than a rectangular array.
18. The method of claim 9 wherein said second hit location being closer to the inside voxel vertex as compared to the outside voxel vertex of the selected o edge.
19. A system for analysis and display of three- dimensional surface structures comprising: first storage means for storing three-dimensional 5 signal patterns representing the value of at least one physical property associated with a three dimensional body at regularly spaced grid locations within said body, a voxel being defined by a set of adjacent in the 0 three-dimensions spaced grid locations, each voxel having voxel edges and voxel faces associated with it; means for defining a plane of pixel positions having preselect position density and 5 corresponding to one plane of said spaced grid locations; means for determining a pixel on said defined pixel plane, said pixel being associated with a voxel that is intersected by a surface structure of interest within said body; processing means for computing coordinate values of the crossection between the surface structure of interest and the intersected voxel edges and for computing the coordinate values of at least one articulating location positioned in the interior of the intersected voxel; second storage means for storing the computed coordinate values along with connectivity information about joining the computed coordinate values into a set of triangles, each triangle having one triangle edge on a voxel face and one triangle vertex at an articulating location in the interior of the intersected voxel; display processor means for receiving data from said second storage means and converting the received data to a display format; and means for interactively displaying surface structures of interest, said display means driven by said display processor means.
20. The system of claim 19 wherein the first storage means receives input from a three-dimensional raw data source.
21. The system of claim 19 wherein said processing means further comprises first generating means for generating the unit normal vectors of each set of triangles associated with a location having computed coordinate values.
22. The system of claim 21 wherein said processing means further comprises second generating means for generating shared normal vectors associated with the voxel locations having computed coordinate values, said shared normal vectors being generated from said triangle normal vectors and said shared normal vectors being stored in said second storage means.
23. The system of claim 22 wherein the display processor means uses said stored shared normal vectors for providing to the display means image data which, when displayed, produces different appearance of different surface structure locations being displayed, said appearance being dependent on said shared normal vectors.
24. The system of claim 19 wherein the processing means further comprises means for identifying which of the voxels sharing a voxel face with the intersected voxel are intersected by the surface structure of interest.
25. The system of claim 24 wherein the processing means further comprises means for generating, for each voxel identified as being intersected, a plurality of index values representing information about the voxel face intersects shared, said plurality of index values being stored in said second storage means.
26. The system of claim 19 wherein the processing means further comprises means for retrieving signal patterns associated with each grid location and means for providing the three-dimensional location associated with each said grid location.
27. The system of claim 19 wherein said processing means comprises means for simultaneously comparing all vertex values of said intersected voxel to a value of said at least one physical property indicative of the surface structure of interest, said comparison used to compute the coordinate values of the crossection independently for each voxel face.
28. The system of claim 19 further comprising variable means for adaptively selecting the position density of the pixels of said defined pixel plane.
29. The system of claim 19 further comprising structure analysis means receiving input from said second storage means and capable of producing on output global structure parameters.
PCT/US1994/004084 1993-04-13 1994-04-13 System and method for surface rendering of internal structures within the interior of a solid object WO1994024640A1 (en)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US046,245 1993-04-13
US08/046,245 US5898793A (en) 1993-04-13 1993-04-13 System and method for surface rendering of internal structures within the interior of a solid object

Publications (1)

Publication Number Publication Date
WO1994024640A1 true WO1994024640A1 (en) 1994-10-27

Family

ID=21942410

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/US1994/004084 WO1994024640A1 (en) 1993-04-13 1994-04-13 System and method for surface rendering of internal structures within the interior of a solid object

Country Status (3)

Country Link
US (1) US5898793A (en)
CA (1) CA2160549A1 (en)
WO (1) WO1994024640A1 (en)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE19603774A1 (en) * 1996-02-02 1997-08-07 Kontron Elektronik Method for determining three-dimensional object properties
EP0794728A1 (en) * 1994-10-27 1997-09-17 Wake Forest University Method and system for producing interactive, three-dimensional renderings of selected body organs having hollow lumens to enable simulated movement through the lumen
WO1999023586A2 (en) * 1997-10-30 1999-05-14 Dr. Baldeweg Gmbh Method and device for processing imaged objects
US5920319A (en) * 1994-10-27 1999-07-06 Wake Forest University Automatic analysis in virtual endoscopy
US5923770A (en) * 1996-09-30 1999-07-13 Siemens Corporate Research, Inc. 3D cardiac motion recovery system using tagged MR images
US6694163B1 (en) 1994-10-27 2004-02-17 Wake Forest University Health Sciences Method and system for producing interactive, three-dimensional renderings of selected body organs having hollow lumens to enable simulated movement through the lumen
CN116804865A (en) * 2023-08-28 2023-09-26 成都飞机工业(集团)有限责任公司 Triaxial automatic programming characteristic identification and tool path generation method

Families Citing this family (64)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH10134208A (en) * 1996-10-31 1998-05-22 Sony Corp Shape data approximating method and plotting device
US8682045B2 (en) * 1997-02-25 2014-03-25 Wake Forest University Health Sciences Virtual endoscopy with improved image segmentation and lesion detection
EP2302594A3 (en) 1998-11-25 2011-04-06 Wake Forest University Virtual endoscopy with improved image segmentation and lesion detection
US6480190B1 (en) * 1999-01-29 2002-11-12 Mitsubishi Electric Research Laboratories, Inc Graphical objects represented as surface elements
US6498607B1 (en) * 1999-01-29 2002-12-24 Mitsubishi Electric Research Laboratories, Inc. Method for generating graphical object represented as surface elements
US6526166B1 (en) * 1999-12-29 2003-02-25 Intel Corporation Using a reference cube for capture of 3D geometry
US6384516B1 (en) 2000-01-21 2002-05-07 Atl Ultrasound, Inc. Hex packed two dimensional ultrasonic transducer arrays
US6772132B1 (en) * 2000-03-02 2004-08-03 Trading Technologies International, Inc. Click based trading with intuitive grid display of market depth
US6754376B1 (en) * 2000-11-22 2004-06-22 General Electric Company Method for automatic segmentation of medical images
US6631202B2 (en) * 2000-12-08 2003-10-07 Landmark Graphics Corporation Method for aligning a lattice of points in response to features in a digital image
US6741246B2 (en) 2001-03-16 2004-05-25 Mitsubishi Electric Research Laboratories, Inc. Hierarchical control point editing of adaptively sampled distance fields
US6724393B2 (en) 2001-03-16 2004-04-20 Mitsubishi Electric Research Labs, Inc. System and method for sculpting digital models
US7034818B2 (en) * 2001-03-16 2006-04-25 Mitsubishi Electric Research Laboratories, Inc. System and method for converting range data to 3D models
US6943789B2 (en) * 2001-03-16 2005-09-13 Mitsubishi Electric Research Labs, Inc Conversion of adaptively sampled distance fields to triangles
US6966837B1 (en) * 2001-05-10 2005-11-22 Best Robert M Linked portable and video game systems
AUPR554301A0 (en) * 2001-06-07 2001-07-12 Howard Florey Institute Of Experimental Physiology And Medicine A map of a property
US6914600B2 (en) * 2001-06-26 2005-07-05 Hewlett-Packard Development Company, L.P. Volumetric warping for voxel coloring on an infinite domain
US20030036083A1 (en) * 2001-07-19 2003-02-20 Jose Tamez-Pena System and method for quantifying tissue structures and their change over time
EP1412909A4 (en) * 2001-07-27 2006-04-19 Virtualscopics Llc System and method for quantitative assessment of joint diseases and the change over time of joint diseases
US20030088177A1 (en) * 2001-09-05 2003-05-08 Virtualscopics, Llc System and method for quantitative assessment of neurological diseases and the change over time of neurological diseases
DE10144004A1 (en) * 2001-09-07 2003-03-27 Philips Corp Intellectual Pty Method of taking three-dimensional measurements for use in measuring physical properties in medical applications, by overlaying three dimensional images
US20030072479A1 (en) * 2001-09-17 2003-04-17 Virtualscopics System and method for quantitative assessment of cancers and their change over time
DE10149556A1 (en) * 2001-10-08 2003-04-24 Siemens Ag Two-dimensional image generation method for medical application, involves defining evaluation surface from three dimensional tomography and selecting surface from predefined data in computer accessible library
US6897863B2 (en) * 2001-11-30 2005-05-24 Caterpillar Inc System and method for hidden object removal
US6867773B2 (en) * 2002-05-13 2005-03-15 Intel Corporation Generating surface normals
US20080262356A1 (en) * 2002-06-07 2008-10-23 Vikram Chalana Systems and methods for ultrasound imaging using an inertial reference unit
US7819806B2 (en) * 2002-06-07 2010-10-26 Verathon Inc. System and method to identify and measure organ wall boundaries
US20090112089A1 (en) * 2007-10-27 2009-04-30 Bill Barnard System and method for measuring bladder wall thickness and presenting a bladder virtual image
US20040127797A1 (en) * 2002-06-07 2004-07-01 Bill Barnard System and method for measuring bladder wall thickness and presenting a bladder virtual image
US20090062644A1 (en) * 2002-06-07 2009-03-05 Mcmorrow Gerald System and method for ultrasound harmonic imaging
US20060025689A1 (en) * 2002-06-07 2006-02-02 Vikram Chalana System and method to measure cardiac ejection fraction
US8221321B2 (en) 2002-06-07 2012-07-17 Verathon Inc. Systems and methods for quantification and classification of fluids in human cavities in ultrasound images
GB2391625A (en) 2002-08-09 2004-02-11 Diagnostic Ultrasound Europ B Instantaneous ultrasonic echo measurement of bladder urine volume with a limited number of ultrasound beams
US7520857B2 (en) * 2002-06-07 2009-04-21 Verathon Inc. 3D ultrasound-based instrument for non-invasive measurement of amniotic fluid volume
US20100036252A1 (en) * 2002-06-07 2010-02-11 Vikram Chalana Ultrasound system and method for measuring bladder wall thickness and mass
US8221322B2 (en) * 2002-06-07 2012-07-17 Verathon Inc. Systems and methods to improve clarity in ultrasound images
JP3760347B2 (en) * 2002-07-30 2006-03-29 株式会社光栄 Program, recording medium, game character drawing method, and game apparatus
US7317456B1 (en) * 2002-12-02 2008-01-08 Ngrain (Canada) Corporation Method and apparatus for transforming point cloud data to volumetric data
US7211046B2 (en) * 2002-12-19 2007-05-01 Ge Medical Systems Global Technology Company Llc Systems and methods for creating reformatted computed tomography images
US7230617B2 (en) * 2003-01-07 2007-06-12 Intel Corporation Voxel center smoothing
US7173615B2 (en) * 2003-01-07 2007-02-06 Intel Corporation Generating a mesh
EP1376381A1 (en) * 2003-02-12 2004-01-02 Agilent Technologies Inc Method and system for data sampling
US20050110793A1 (en) * 2003-11-21 2005-05-26 Steen Erik N. Methods and systems for graphics processing in a medical imaging system
EP1577837A1 (en) * 2004-03-17 2005-09-21 Deutsches Krebsforschungszentrum Stiftung des öffentlichen Rechts 3D cone beam reconstruction
US7633502B2 (en) * 2004-05-19 2009-12-15 Boston Scientific Scimed, Inc. System and method for graphically representing anatomical orifices and vessels
US7571081B2 (en) * 2004-07-15 2009-08-04 Harris Corporation System and method for efficient visualization and comparison of LADAR point data to detailed CAD models of targets
JP4783100B2 (en) * 2005-09-12 2011-09-28 独立行政法人理化学研究所 Method of converting boundary data into in-cell shape data and its conversion program
JP4855853B2 (en) * 2006-07-05 2012-01-18 富士通株式会社 Analysis device, computer control method, and model creation program
US8217937B2 (en) * 2007-03-28 2012-07-10 The Aerospace Corporation Isosurfacial three-dimensional imaging system and method
US8167803B2 (en) * 2007-05-16 2012-05-01 Verathon Inc. System and method for bladder detection using harmonic imaging
US8133181B2 (en) * 2007-05-16 2012-03-13 Verathon Inc. Device, system and method to measure abdominal aortic aneurysm diameter
US8131076B2 (en) * 2007-06-14 2012-03-06 Siemens Aktiengesellschaft Editing of pre-segmented images using seeds derived from contours
US8243334B2 (en) 2008-06-06 2012-08-14 Virginia Venture Industries, Llc Methods and apparatuses for printing three dimensional images
DE102008032006B4 (en) * 2008-07-07 2017-01-05 Siemens Healthcare Gmbh Method for controlling the image recording in an image recording device, and an image recording device
US8225998B2 (en) * 2008-07-11 2012-07-24 Es&S Innovations Llc Secure ballot box
US8217939B1 (en) * 2008-10-17 2012-07-10 Ngrain (Canada) Corporation Method and system for calculating visually improved edge voxel normals when converting polygon data to voxel data
EP3427687A1 (en) * 2009-05-14 2019-01-16 Covidien LP Automatic registration technique
US8655033B2 (en) * 2009-10-28 2014-02-18 General Electric Company Iterative reconstruction
US9082218B2 (en) * 2012-03-23 2015-07-14 Roblox Corporation Method for manipulating three-dimensional voxel data for on-screen visual
CN104658034B (en) * 2013-11-18 2019-03-01 清华大学 The fusion method for drafting of CT image data
US11006852B2 (en) 2017-12-11 2021-05-18 Covidien Lp Systems, methods, and computer-readable media of estimating thoracic cavity movement during respiration
US11620478B2 (en) * 2019-02-06 2023-04-04 Texas Instruments Incorporated Semantic occupancy grid management in ADAS/autonomous driving
USD1011573S1 (en) 2021-03-18 2024-01-16 Milwaukee Electric Tool Corporation Lighting apparatus
US11810255B2 (en) * 2021-05-28 2023-11-07 Adobe Inc. Swept volume determination techniques

Citations (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4525858A (en) * 1983-01-03 1985-06-25 General Electric Company Method and apparatus for reconstruction of three-dimensional surfaces from interference fringes
US4710876A (en) * 1985-06-05 1987-12-01 General Electric Company System and method for the display of surface structures contained within the interior region of a solid body
US4719585A (en) * 1985-08-28 1988-01-12 General Electric Company Dividing cubes system and method for the display of surface structures contained within the interior region of a solid body
US4729098A (en) * 1985-06-05 1988-03-01 General Electric Company System and method employing nonlinear interpolation for the display of surface structures contained within the interior region of a solid body
US4751643A (en) * 1986-08-04 1988-06-14 General Electric Company Method and apparatus for determining connected substructures within a body
US4791567A (en) * 1986-09-15 1988-12-13 General Electric Company Three dimensional connectivity system employing an equivalence schema for determining connected substructures within a body
US4821210A (en) * 1987-04-02 1989-04-11 General Electric Co. Fast display of three-dimensional images
US4821213A (en) * 1986-12-19 1989-04-11 General Electric Co. System for the simultaneous display of two or more internal surfaces within a solid object
US4831528A (en) * 1987-11-09 1989-05-16 General Electric Company Apparatus and method for improvement of 3D images derived from tomographic data
US4879668A (en) * 1986-12-19 1989-11-07 General Electric Company Method of displaying internal surfaces of three-dimensional medical images
US4885688A (en) * 1987-11-25 1989-12-05 General Electric Company Minimization of directed points generated in three-dimensional dividing cubes method
US4984157A (en) * 1988-09-21 1991-01-08 General Electric Company System and method for displaying oblique planar cross sections of a solid body using tri-linear interpolation to determine pixel position dataes
US5113357A (en) * 1989-05-18 1992-05-12 Sun Microsystems, Inc. Method and apparatus for rendering of geometric volumes
US5166876A (en) * 1991-01-16 1992-11-24 General Electric Company System and method for detecting internal structures contained within the interior region of a solid object

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0463700B1 (en) * 1990-06-29 1997-09-03 Philips Electronics Uk Limited Method of and apparatus for generating an image

Patent Citations (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4525858A (en) * 1983-01-03 1985-06-25 General Electric Company Method and apparatus for reconstruction of three-dimensional surfaces from interference fringes
US4710876A (en) * 1985-06-05 1987-12-01 General Electric Company System and method for the display of surface structures contained within the interior region of a solid body
US4729098A (en) * 1985-06-05 1988-03-01 General Electric Company System and method employing nonlinear interpolation for the display of surface structures contained within the interior region of a solid body
US4719585A (en) * 1985-08-28 1988-01-12 General Electric Company Dividing cubes system and method for the display of surface structures contained within the interior region of a solid body
US4751643A (en) * 1986-08-04 1988-06-14 General Electric Company Method and apparatus for determining connected substructures within a body
US4791567A (en) * 1986-09-15 1988-12-13 General Electric Company Three dimensional connectivity system employing an equivalence schema for determining connected substructures within a body
US4879668A (en) * 1986-12-19 1989-11-07 General Electric Company Method of displaying internal surfaces of three-dimensional medical images
US4821213A (en) * 1986-12-19 1989-04-11 General Electric Co. System for the simultaneous display of two or more internal surfaces within a solid object
US4821210A (en) * 1987-04-02 1989-04-11 General Electric Co. Fast display of three-dimensional images
US4831528A (en) * 1987-11-09 1989-05-16 General Electric Company Apparatus and method for improvement of 3D images derived from tomographic data
US4885688A (en) * 1987-11-25 1989-12-05 General Electric Company Minimization of directed points generated in three-dimensional dividing cubes method
US4984157A (en) * 1988-09-21 1991-01-08 General Electric Company System and method for displaying oblique planar cross sections of a solid body using tri-linear interpolation to determine pixel position dataes
US5113357A (en) * 1989-05-18 1992-05-12 Sun Microsystems, Inc. Method and apparatus for rendering of geometric volumes
US5166876A (en) * 1991-01-16 1992-11-24 General Electric Company System and method for detecting internal structures contained within the interior region of a solid object

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
NIELSON et al., "The Asymptotic Decider: Resolving the Ambiguity in Matching Cubes", COMPUTER SCIENCE, ARIZONA STATE UNIVERSITY, 1991, pp. 83-91. *

Cited By (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6366800B1 (en) 1994-10-27 2002-04-02 Wake Forest University Automatic analysis in virtual endoscopy
US7792565B2 (en) 1994-10-27 2010-09-07 Wake Forest University Health Sciences Method and system for producing interactive, three-dimensional renderings of selected body organs having hollow lumens to enable simulated movement through the lumen
US6694163B1 (en) 1994-10-27 2004-02-17 Wake Forest University Health Sciences Method and system for producing interactive, three-dimensional renderings of selected body organs having hollow lumens to enable simulated movement through the lumen
US8275446B2 (en) 1994-10-27 2012-09-25 Wake Forest University Health Sciences Automatic analysis in virtual endoscopy
US5920319A (en) * 1994-10-27 1999-07-06 Wake Forest University Automatic analysis in virtual endoscopy
US7853310B2 (en) 1994-10-27 2010-12-14 Wake Forest University Health Sciences Automatic analysis in virtual endoscopy
EP0794728A1 (en) * 1994-10-27 1997-09-17 Wake Forest University Method and system for producing interactive, three-dimensional renderings of selected body organs having hollow lumens to enable simulated movement through the lumen
US6272366B1 (en) 1994-10-27 2001-08-07 Wake Forest University Method and system for producing interactive three-dimensional renderings of selected body organs having hollow lumens to enable simulated movement through the lumen
US7149564B2 (en) 1994-10-27 2006-12-12 Wake Forest University Health Sciences Automatic analysis in virtual endoscopy
EP0794728A4 (en) * 1994-10-27 1998-03-11 Univ Wake Forest Method and system for producing interactive, three-dimensional renderings of selected body organs having hollow lumens to enable simulated movement through the lumen
US6909913B2 (en) 1994-10-27 2005-06-21 Wake Forest University Health Sciences Method and system for producing interactive three-dimensional renderings of selected body organs having hollow lumens to enable simulated movement through the lumen
DE19603774A1 (en) * 1996-02-02 1997-08-07 Kontron Elektronik Method for determining three-dimensional object properties
US5923770A (en) * 1996-09-30 1999-07-13 Siemens Corporate Research, Inc. 3D cardiac motion recovery system using tagged MR images
WO1999023586A2 (en) * 1997-10-30 1999-05-14 Dr. Baldeweg Gmbh Method and device for processing imaged objects
WO1999023586A3 (en) * 1997-10-30 1999-07-08 Baldeweg Gmbh Dr Method and device for processing imaged objects
US6734847B1 (en) 1997-10-30 2004-05-11 Dr. Baldeweg Gmbh Method and device for processing imaged objects
CN116804865B (en) * 2023-08-28 2023-12-08 成都飞机工业(集团)有限责任公司 Triaxial automatic programming characteristic identification and tool path generation method
CN116804865A (en) * 2023-08-28 2023-09-26 成都飞机工业(集团)有限责任公司 Triaxial automatic programming characteristic identification and tool path generation method

Also Published As

Publication number Publication date
CA2160549A1 (en) 1994-10-27
US5898793A (en) 1999-04-27

Similar Documents

Publication Publication Date Title
US5898793A (en) System and method for surface rendering of internal structures within the interior of a solid object
US5166876A (en) System and method for detecting internal structures contained within the interior region of a solid object
US7920734B2 (en) Method for visualization of plaque deposits from 3D image data records of vessel structures
JP2744490B2 (en) Apparatus and method for displaying two-dimensional image of object internal structure surface
US6990231B2 (en) Method and apparatus for forming and displaying projection image from a plurality of sectional images
Sherbondy et al. Exploring connectivity of the brain's white matter with dynamic queries
US5412563A (en) Gradient image segmentation method
CA2278672C (en) Automatic analysis in virtual endoscopy
US4719585A (en) Dividing cubes system and method for the display of surface structures contained within the interior region of a solid body
EP1365356B1 (en) Semi-automatic segmentation algorithm for pet oncology images
JP2885842B2 (en) Apparatus and method for displaying a cut plane in a solid interior region
US7576740B2 (en) Method of volume visualization
EP0204225A2 (en) System and method for the display of surface structures contained within the interior region of a solid body
US4953087A (en) Three-dimensional images obtained from tomographic data having unequally spaced slices
Mankovich et al. Three-dimensional image display in medicine
Udupa 3D imaging: principles and approaches
WO2004095378A1 (en) Combined 3d and 2d views
EP0373854B1 (en) Apparatus and method for detecting internal structures contained within the interior region of a solid object
Mullick et al. Confocal volume rendering: fast, segmentation-free visualization of internal structures
JP2005525863A (en) Medical inspection system and image processing for integrated visualization of medical data
Rusinek et al. Quantitative and qualitative comparison of volumetric and surface rendering techniques
EP0501812A2 (en) Simulated surgery using display list surface data
de Araujo Buck et al. 3-D segmentation of medical structures by integration of raycasting with anatomic knowledge
Barrett Interactive measurement and display of three-and four-dimensional anatomy
Wood et al. Artifacts and illusions in surface and volume rendering

Legal Events

Date Code Title Description
AK Designated states

Kind code of ref document: A1

Designated state(s): CA

AL Designated countries for regional patents

Kind code of ref document: A1

Designated state(s): AT BE CH DE DK ES FR GB GR IE IT LU MC NL PT SE

DFPE Request for preliminary examination filed prior to expiration of 19th month from priority date (pct application filed before 20040101)
121 Ep: the epo has been informed by wipo that ep was designated in this application
WWE Wipo information: entry into national phase

Ref document number: 2160549

Country of ref document: CA

122 Ep: pct application non-entry in european phase