Method for modelling a surface and device for implementing same

Abstract

A method for obtaining a model of a surface including the steps of obtaining measurements of geometrical data concerning specific points on the surface, making a grid of the surface, with the grid passing through said points, memorizing, at an address which is specific to each node of the grid, the coordinates of the node, the number of satellites of the node, information for access to the addresses of said satellites and thereafter to information which relates to them, and geometrical data which may be associated with said node; for each node, defining a local roughness index obtained from a weighted sum of the current coordinates of the node and its satellites, defining the sum of an overall roughness index representing the sum of all the local roughness indices, and of an overall index of the infringement of said geometrical data, iteratively adjusting the coordinates of indefinite nodes, by using at each adjustment the sum of a weighted combination of current node neighbor coordinates and of a combination of geometrical data associated with the node, in order to minimize said sum, and creating a model of the surface on the basis of the adjusted coordinates.

Description

SUMMARY OF THE INVENTION

The present invention is generally concerned with investigations within a three-dimensional body and is particularly concerned with a new process of obtaining a representation of a surface within a body on the basis of a limited set of known geometrical data relating to said surface.

Investigations within three-dimensional bodies are of major concern in geology, geophysics and biology.

In geophysics, for example, it is necessary to obtain as accurate as possible representations of surfaces situated, for example, at the interface between two areas of different kinds or with different properties, on the basis of data obtained from prospecting, or during exploitation of an underground resource.

How effectively the three-dimensional body is investigated, especially in prospecting for oil or other underground resources, depends on the accuracy with which this type of surface can be reconstituted and represented.

In medicine and biology there are various known processes for obtaining representations of cross-sections through a living body or the like and the aim is also to reconstitute and represent in three dimensions a surface situated at the interface between two media, for example the contours of a bodily organ.

There already exist a number of modeling techniques for obtaining a representation of a surface within a three-dimensional body that is to be exploited or treated. Particularly worthy of mention are the Bezier interpolation method and the spline functions method (for more details reference should be had to Geometric Modeling, M. E. Mortenson, Ed. John Wiley, 1985).

However, all these known techniques are ill suited to handling heterogeneous data, by which is meant data relating to the coordinates of points on the surface, to fixed orientations of the surface, to geometrical links between two separate points, etc. Also, these techniques are ill-suited to a wide spread of surfaces or to surface anomalies such as folds and breaks, and to discontinuous surfaces. To be more precise, in some cases the functions used in these techniques are not convergent or have no solution or have no single solution.

Finally, the known techniques are not able to take into account the concept of the degree of certainty with which some surface data is known.

The present invention is directed to alleviating the drawbacks of the prior art processes and to proposing a process that can allow for non-homogeneous and/or highly diverse geometrical data and for the anomalies mentioned above and other anomalies, providing in each case a single and unambiguous model of the surface.

Another object of the invention is to propose a process that can use geometrical data and at the same time data as to the degree of certainty or accuracy of said data.

To this end, the present invention firstly consists in a process for modeling a surface representing for example the interface between two areas of different kinds or with different properties in a three-dimensional body such as a geological formation or a living body, characterized in that it comprises the steps of:

obtaining by means of measuring apparatus a set of geometrical data relating to the surface and associated with respective points on said surface; meshing the surface so that all said points are a subset of the nodes of the mesh; storing at a specific memory address for each node of the mesh, the following data: the coordinates of the node in question, the number of satellite nodes of the node in question, data providing access to the specific addresses of said satellite nodes and consequently to the data relating thereto, if necessary, geometrical data associated with said node in question, for each node of the mesh, defining a local roughness index derived from a weighted sum of the actual coordinates of the node and of its satellites, defining the sum of a global roughness index obtained by summing the local roughness indices associated with each node and a global index of violation of said geometrical data, fitting the coordinates of each node for which the precise coordinates are not known by an iterative method in which for each step of the iteration there are added a weighted combination of the actual coordinates of the satellite and of the satellites of the satellites of said node and a combination of the geometrical data associated with said node, in such a way as to minimize said sum, and creating a representation of the surface from the fitted coordinates of each node.

The invention also concerns a device for modeling a surface representing for example the interface between two areas of different kinds or with different properties in a three-dimensional body such as a geological formation or a living body, characterized in that it comprises:

means for obtaining by means of measuring apparatus a set of geometrical data relating to the surface and associated with respective points on said surface; means for meshing the surface so that all said points are a subset of the nodes of the mesh; means for storing at a specific memory address for each node of the mesh, the following data: the coordinates of the node in question, the number of satellite nodes of the node in question, data providing access to the specific addresses of said satellite nodes and consequently to the data relating thereto, if necessary, geometrical data associated with said node in question, calculation means for fitting the coordinates of each node for which the precise coordinates are not known by an iterative method in which for each step of the iteration there are added a weighted combination of the actual coordinates of the satellite and of the satellites of the satellites of said node and a combination of the geometrical data associated with said node, in such a way as to minimize a sum of a global roughness index obtained by summing local roughness indices associated with each node and each derived from a weighted sum of the actual coordinates of the node and of its satellites and a global index of violation of said geometrical data, means for creating a representation of the surface from the fitted coordinates of each node.

These and other objects of the present invention will become more readily apparent from the detailed description given hereinafter. However, it should be understood that the detailed description and specific examples, while indicating preferred embodiments of the invention, are given by way of illustration only, since various changes and modifications within the spirit and scope of the invention will become apparent to those skilled in the art from this detailed description.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention will emerge more clearly from the following detailed description of one preferred embodiment of the invention given by way of non-limiting example and with reference to the appended drawings, in which:

FIG. 1 shows a typical graph incorporating meshing used to model a surface, this graph being defined by the set of all the sides of facets whose vertices constitute the set .OMEGA. described later, and also shows arrows symbolically representing displacements of nodes of the mesh applied interactively by the user to shape the surface;

FIGS. 2(a) through 2(d) show four instances of controlling the shape of the surface using four different types of vector constraint; the vectors .DELTA..sub..lambda..mu. are assumed to be known directly in FIGS. 2(a) through 2(c) and known because they are orthogonal to a given vector V in FIG. 2(d);

FIG. 3 shows a typical complex geological surface (salt dome) modeled using the process of the present invention;

FIG. 4 shows a typical complex biological surface (embryo brain) modeled using the process of the present invention from discrete data obtained from a succession of cross-sections;

FIG. 5 shows the organization of the memory associated with an atom A(k) of the kernel .phi..sub.k ;

FIG. 6 shows a first example of memorizing an orbit of a given atom in an array with n.sub.k items;

FIG. 7 shows a second example of memorizing an orbit of a given atom in n.sub.k arrays each of two items;

FIG. 8 shows the memorization of a geometrical datum or constraint;

FIG. 9 shows the method of memorizing a set of geometrical data or constraints relating to a given atom A(k);

FIG. 10 shows the memorization in an array of all of the atoms relating to a surface S; and

FIG. 11 shows the chaining of access to said set of atoms;

FIG. 12 shows an apparatus embodiment of the present invention.

Identical or similar elements or parts are designated by the same reference numbers in all the figures.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENT

Reference will be made later to a number of publications whose respective contents are deemed to be included by way of reference into this description.

(1) A method of Bivariate interpolation and smooth surface fitting based on local procedures, H. AKIMA, J. ACM 17.1, 1974;

(2) Shape reconstruction from planar cross-sections, J. D. BOISSONNAT, ACM Transactions on Graphics, vol. No. 2, 1986;

(3) Machine contouring using minimum curvature, I. C. BRIGGS, Geophysics 17.1, 1974;

(4) Triangular Bernstein-Bezier patches, G. FARIN, Computer-aided geometric design, 1986;

(5) Primitives for the manipulation of general subdivisions and the computation of Voronoi diagrams, L. GUIBAS and J. STOLFI, ACM transactions on Graphis, vol. No. 2, 1985;

(6) Three-dimensional graphic display of disconnected bodies, J. L. MALLET, Mathematical Geology, vol. 20, no. 8, 1988;

(7) Geometrical modeling and geostatistics, J. L. MALLET, Proceedings of the Third International Geostatistics Congress, Kluwer Academic Publishers, 1989;

(8) Discrete Smooth Interpolation, J. L. MALLET, ACM Transactions on Graphics, April issue, 1989;

(9) Geometric Modeling, M. E. MORTENSON, John Wiley, 1985.

Also, the following description of one preferred embodiment of the invention is followed by section 2, infra explaining a method of reconstituting a surface derived from the method explained in reference document (8) and on which the process of the present invention is based.

1. DESCRIPTION OF A PROCESS IN ACCORDANCE WITH THE INVENTION FOR ENCODING AND FITTING COMPLEX SURFACES

1.1 Introduction

1.1.1 Nature of the process

An industrial process described hereafter is used to encode and fit complex surfaces encountered, for example, when modeling:

interfaces between different domains, interfaces between geological layers, interfaces between biological organs, etc.

The result of this encoding phase may be converted into pictures and/or solid models (molded from plastics material, for example) by a graphic computer or any other device. In geology, for instance, these pictures and/or solid models are used to optimize the prospecting, exploitation and management of earth resources such as:

oil deposits, mineral deposits, water deposits, etc.

The process described hereafter consists of two complementary sub-processes:

1. a sub-process to encode a surface, 2. a sub-process to fit a surface to the data currently available. In geology, for example, these data may be: the exact or approximate location of points in the intersection of a given surface and a plane, the exact or approximate location of the intersection of a well and a geological surface corresponding to the interface between two geological layers, the exact or approximate location of the plane tangent to a geological surface corresponding to the interface between two geological layers, the exact or approximate throw vector of a fault breaking the geological surface corresponding to the interface between two geological layers, etc.

The encoding process is designed to optimize the fitting process based on a new version of the DSI method (cf. [8]) especially developed for the geometrical modeling of surfaces. This new version of the DSI method, which has not yet been published, is presented in section 2, infra.

1.1.2 Dividing a surface into facets

Let S be the surface to model. As shown in FIG. 1 and as explained in [6] and [7], this surface will be approximated by a set of triangular and/or polygonal facets. These facets are characterized by:

The coordinates of their vertices. These vertices constitute a finite set of N points {.phi..sub.1, . . . ,.phi.N } regularly distributed on S and numbered from 1 to N. The coordinates of the k.sup.th vertex .phi..sub.k will be denoted (.phi..sub.k.sup.x,.phi..sub.k.sup.y,.phi..sub.k.sup.z). Their edges constituting the links between pairs of vertices (.phi..sub.i,.phi..sub.j).

The set of edges of the triangles and/or polygons composing the surface S constitute a graph G whose nodes are the vertices of the triangles and/or polygons. For this reason and for the sake of simplicity, "node" is used hereafter as a synonym of "vertex of triangles and/or polygons".

node.uparw.{vertex of triangles and/or polygons} 1.1.3 Preliminary definitions

Set of satellites .LAMBDA.(k)

Let .phi..sub.k be the k.sup.th node of a surface S. The "set of satellites" of .phi..sub.k is the set .LAMBDA.(k) Of the nodes .phi..sub.j different from .phi..sub.k and such that (.phi..sub.j,.phi..sub.k) iS an edge of triangle or polygon. In other words, the set .LAMBDA.(k) is such that: ##EQU1##

As shown above, the number of satellites attached to the node .phi..sub.k is denoted n.sub.k.

The orbit concept discussed in the articles [6] and [7] has not yet been subjected to any particular encoding process; an important feature of the invention is its proposed method of encoding .LAMBDA.(k).

Neighborhood N(k)

Let .phi..sub.k be the k.sup.th node of a surface S. The neighborhood N(k) of .phi..sub.k is the subset N(k) of nodes equal to the union of .phi..sub.k and its satellites. In other words:

N(k)={.phi..sub.k }.orgate..LAMBDA.(k)

A neighbor of node .phi..sub.k is any node belonging to the set N(k). Thus:

.phi..sub.j .epsilon.N(k) {.phi..sub.j =Neighbor of .phi..sub.k }

The neighborhood concept discussed in the articles [7] and [8] has not yet been subjected to any particular encoding process; an important feature of the invention is its proposed method of encoding N(k).

Atom A(R) and associated kernel

Let .phi..sub.k be the k.sup.th node of a surface S. The atom of kernel .phi..sub.k A(R) is the pair (.phi..sub.k,.LAMBDA.*(k)) composed of:

the node .phi..sub.k the set .LAMBDA.*(K) of the addresses or the codes {A*(j1), . . . ,A*(jn.sub.k)} allowing the set of the corresponding atoms {A(j1), . . . ,A(jn.sub.k)} whose kernels are the satellites of .phi..sub.k to be retrieved from a memory.

In other words: ##EQU2##

The atom concept discussed in the article [6] has not yet been subjected to any particular encoding process; an important feature of the invention is its proposed method of encoding A(k).

Constraints and types of constraints

A constraint attached to an atom A(k) is any condition concerning the location of the node .phi..sub.k corresponding to the kernel of this atom. For example, in geology, the following constraints are of prime importance to fitting a surface to precise or imprecise data observed or measured by a manual or geophysical technique:

The coordinates (.phi..sub.k.sup.x,.phi..sub.k.sup.y,.phi..sub.k.sup.z) of the node .phi..sub.k are known exactly. Such a constraint is called a control node type constraint and may be used in geology to encode the intersection of the modeled surface with a well. The coordinates (.phi..sub.k.sup.x,.phi..sub.k.sup.y,.phi..sub.k.sup.z) of the node .phi..sub.k are not known exactly but with a given certainty factor. Such a constraint is called a fuzzy control node type constraint and may be used to encode the approximate location of a node of the modeled surface measured approximately by a given process. For example, seismic data in geology, and data corresponding to sections of organs observed with a microscope or using ultrasound techniques or with a scanner in biology are of this type. The vector (.phi..sub..lambda. -.phi..sub.k) joining the node .phi..sub.k to another node .phi..sub..lambda. is approximately known with a given certainty factor. Such a constraint is called a fuzzy vectorial link type constraint and is used in geology, for example, to encode the throw vector of a fault, the nodes .phi..sub.k and .phi..sub..lambda. being on opposite sides of the fault (see FIG. 2). The vector V orthogonal to the tangent plane to the surface S at node .phi..sub.k is approximately known with a given certainty factor. Such a constraint is called a fuzzy vector normal type constraint and is used in geology, for example, to encode the plane tangent to layer interfaces measured by a manual or geophysical technique (sounding). Etc.

The certainty factors mentioned above are positive numbers whose values are assumed to be proportional to the a priori confidence that can be attached to the corresponding data. Accounting for these fuzzy constraints constitutes one of the characteristics of the invention and is described in detail in section 2 (see sections 2.5 and 2.7.6), infra).

Notation

To simplify, no distinction will be drawn between the number of any node of a surface and the node itself. This implies that:

.LAMBDA.(k)=set of satellites {.phi..sub.j1, . . . ,.phi..sub.jn.sbsb.k } of .phi..sub.k or set of corresponding numbers {j1, . . . ,jn.sub.k }, N(k)=set of neighbors {.phi..sub.k,.phi..sub.j1, . . . ,.phi..sub.jn.sbsb.k } of .phi..sub.k or set of corresponding numbers {k,j1, . . . ,jn.sub.k }. 1.1.4 State of the art

Classical methods

A complete description of the "state of the art" concerning the modeling of complex surfaces is given in reference [7]. This article explains precisely why the classical methods used in automatic mapping and in Computer Aided Design (CAD) are not relevant to the modeling of complex surfaces, as might otherwise be expected. The basics of the DSI method (also known as the DSI/DSA method) which overcomes these difficulties are briefly introduced at the end of reference [7].

The DSI method

This method is based on the notion of "roughness" of a surface S composed of triangular and/or polygonal facets. This notion of "roughness" R(.phi..vertline.k) at node .phi..sub.k of a surface S is defined as follows (See ref. [7], [8] and appendix): ##EQU3##

The positive or null coefficients {v.sup..alpha. (k)} occurring in this definition are weighting coefficients chosen by the user. Among the infinity of possible choices, one of the most interesting, called "harmonic weighting", consists of choosing: ##EQU4##

The local roughnesses so defined are then combined in a global roughness R(.phi.): ##EQU5##

The positive or null coefficients .mu.(k) occurring in this definition are used to modulate the contribution of the local roughness (R(.phi..vertline.k) to the global roughness R(.phi.). If a uniform weighting is used, for example:

.mu.(k)=1 k

Among other possible choices for the coefficients .mu.(k), the following weightings can be used, where m is a given positive constant: ##EQU6##

The principle of the DSI method consists in computing the location of each node .phi..sub.k of the surface S to model in order to render it as smooth as possible while respecting the data and constraints governing the shape of the surface. The paper [8] is devoted to a mathematical presentation of the DSI method and proposes a technique for minimizing the global roughness while respecting the data and constraints.

Interactive version Of the DSI method

The method proposed in [8] for minimizing the global roughness is not well suited to interactive use on an electronic computer (or microprocessor), which is why section 2, infra proposes a new method on which the invention is based. As mentioned previously, the object of the invention is to propose an encoding and fitting process for the surface S based on this new version of the DSI method.

1.2 Examples of images of objects encoded and fitted with the process

1.2.1 Geological example

Oil prospecting is generally performed indirectly in the sense that one is looking for geological layers which could have acted as a trap for oil; for example, deformations of salt layers called "salt domes" are generally excellent oil traps. In order to locate the oil precisely, it is necessary to know the exact shape of the "trap" layers and unfortunately, in many cases, until now there was no practical process usable to model it efficiently. The salt dome case is generally regarded as the most complex and that is why to prove the efficiency of this encoding and fitting process FIG. 3 shows an image of such a surface obtained with it.

1.2.2 Biological example

Medical data acquisition processes (scanner or ultrasound) are mathematically identical to seismic techniques used in oil prospecting. For example, the ultrasound technique is identical to seismic reflection and scanners are based on tomography as in seismic tomography. This is one reason why this encoding and fitting process may be used to model the skin of biological organs. FIG. 4 shows a model of the brain of a 7 millimeters long human embryo obtained using this encoding and fitting method. Note that until now organs this small (the brain is 1 millimeter long) could only be seen in the form of microscopic sections.

1.3 Encoding process for a surface S

1.3.1 Memory and memory address concepts

The term "memory" denotes any electronic computing device able to record data, codes or partial results and to retrieve them on request for use in subsequent operations. Moreover, the "address" of a memory is any code identifying the physical location of the memory inside the computer.

1.3.2 Encoding an atom A(k)

The following process or one of its variants may be used to encode and store each atom A(k) in a set of memories:

Encoding process

As suggested in FIG. 5, a contiguous memory area beginning at address A*(k) stores consecutively the following data relating to the atom A(k) corresponding to the node .phi..sub.k :

1. Kernel= either the coordinates (.phi..sub.k.sup.x,.phi..sub.k.sup.y,.phi..sub.k.sup.z) of the node .phi..sub.k, or the address or the code or any other data allowing the coordinates (.phi..sub.k.sup.x,.phi..sub.k.sup.y,.phi..sub.k.sup.z) of the node .phi..sub.k to be retrieved. 2. nb.sub.-- sat= either the number n.sub.k of satellites linked to the node .phi..sub.k, or the address or the code or any other data allowing the number n.sub.k of satellites linked to the node .phi..sub.k to be retrieved. 3. Sat=

Address of a memory area containing the addresses or the access codes {A*(j1), . . . ,A*(jn.sub.k)} allowing the n.sub.k atoms {A(j1), . . . ,A(jn.sub.k)} whose kernels {.phi..sub.j1, . . . ,.phi..sub.jn.sbsb.k } are the satellites of the node .phi..sub.k to be retrieved.

4. Const=

Encoded constraints (see section 1.1.3) attached to the node .phi..sub.k. This encoding will be explained in detail in section 1.1.3.

5. Info=

Memory area of variable size that can contain complementary data concerning the atom A(k). This data is specific to each particular application; for instance, in oil prospecting, if A(k) is an atom of the surface separating two adjacent layers C1 and C2, then it may be interesting to store in the Info field:

a list of physical properties (porosities, seismic velocities, etc) of the layer C1 at node .phi..sub.k corresponding to the atom A(k), a list of physical properties (porosities, seismic velocities, etc) of the layer C2 at node .phi..sub.k corresponding to the atom A(k), a list of geological properties (geological facies, presence of oil, etc) of the layer C1 at node .phi..sub.k corresponding to the atom A(k), a list of geological properties (geological facies, presence of oil, etc) of the layer C2 at node .phi..sub.k corresponding to the atom A(k), etc.

Variant 1

The memory area in which are stored the addresses or the access codes {A*(j1), . . . ,A*(jn.sub.k)} for retrieving the n.sub.k atoms {A(j1), . . . ,A(jn.sub.k)} whose kernels are the satellites {.phi..sub.j1, . . . ,.phi..sub.jn.sbsb.k } of the node .phi..sub.k may be structured in two ways:

1. As shown in FIG. 6, the first way consists in using an array composed of n.sub.k consecutive memories containing the addresses {A*(j1), . . . ,A*(jn.sub.k)}. 2. The second way (FIG. 7) consists in using n.sub.k arrays each containing two consecutive memories such that for the .alpha..sup.th array: the first memory contains the address A*(j.sub..alpha.) of the atom corresponding to the .alpha..sup.th satellite of the atom A(k), the second memory contains the address of the (.alpha.+1).sup.th array. If .alpha. is equal to the number n.sub.k of satellites (nb.sub.-- sat=n.sub.k), then this second memory is either unused or contains a code specifying that the last satellite of A(k) has been reached.

Variant 2

It is possible not to store the nb.sub.-- sat field but in this case, if the list of addresses or access codes {A*(j1), . . . ,A*(jn.sub.k)} iS coded as a single array (see variant 1), it is necessary to add a memory area A*(n.sub.k +1) to store a code indicating that the list .LAMBDA.*(k) is finished.

Variant 3

The Sat field has to be made big enough to contain the addresses or access codes {A*(j1), . . . ,A*(jn.sub.k)}.

Variant 4

Many variants are possible depending on the size of the Info field and on how it is partitioned. However, it is obvious that all these variants have no influence on the use of the atom A(k) in connection with the fitting algorithm based on the DSI method described in section 2, infra.

Variant 5

Other fields can be added to the atom A(k) described above. However these other fields have no influence on the use of the atom A(k) in connection with the fitting algorithm based on the DSI method described in section, infra.

Variant 6

For any variant, what is really important for an efficient implementation of the DSI method is that the address or the access code A*(k) of a given atom A(k) should be the most direct path to the atoms whose kernels are the satellites of the kernel of A(k).

The main object of the method of encoding A(k) in a memory described above is to allow direct access to the satellites.

1.3.3 Encoding the Constraints

Each constraint relative to an atom A(k) will be stored in a specific memory area whose description follows:

Process for encoding a constraint

As shown in FIG. 8, the encoding of a particular constraint attached to a given atom is performed in a memory area made up of two fields, Constraint and Next. These two fields are used as follows:

Constraint=

Memory area containing:

either the data related to the encoded constraint (see below), or the address or the access code of a memory area containing the data related to the encoded constraint.

The data related to the encoded constraint must include a code identifying the type of the encoded constraint (see section 1.1.3).

Next=

Memory area containing:

either a code indicating that the current constraint is the last constraint attached to the atom A(k), or the address of a memory area containing the next constraint (see below).

Process for storing all constraints attached to an atom

As described in section 1.3.2, encoding an atom A(k) requires a Const field that will be used to store:

either a code specifying that the atom A(k) has no constraint, or the address of a memory area containing the first constraint related to the atom A(k).

If the address of the n.sup.th constraint related to the atom A(k) is stored in the Next field of the (n-1).sup.th constraint (n>1), then, as shown in FIG. 9, this implies that all the constraints related to A(k) are attached to A(k). This encoding process allows new constraints to be added or old constraints to be removed without modifying the organization of the other data.

Examples of data relating to a given constraint

Section 1.1.3 gives three examples of constraint types which are of prime importance for the encoding and fitting of complex surfaces encountered in biology and geology. For each of these types of constraints data has to be attached to the Constraint field previously described and shown in FIG. 9:

For a fuzzy control node type constraint, in addition to the constraint type, it is necessary to store at least the following data: the approximate coordinates (.phi..sub.k.sup.x,.phi..sub.k.sup.x,.phi..sub.k.sup.z) of the node .phi..sub.k corresponding to the kernel of the atom A(k) to which the constraint is attached, the certainty factor (positive number) proportional to the confidence that can be placed in the approximate coordinates (.phi..sub.k.sup.x,.phi..sub.k.sup.y,.phi..sub.k.sup.z). For a fuzzy vectorial link type constraint, in addition to the type, it is necessary to store at least the following data: the address of the node .phi..sub..lambda. linked to the kernel .phi..sub.k of the atom A(k) to which the constraint is attached, the three approximate coordinates of the vector (.phi..sub..lambda. -.phi..sub.k) joining the node .phi..sub.k to the node .phi..sub..lambda., the certainty factor (positive number) proportional to the confidence that can be placed in the approximate coordinates of the vector (.phi..sub..lambda. -.phi..sub.k). For a fuzzy vector normal type constraint, in addition to the type, it is necessary to store at least the following data: the three approximate components of the vector V orthogonal to the surface S at the node .phi..sub.k corresponding to the kernel of the atom A(k) to which the constraint is attached, the certainty factor (positive number) proportional to the confidence that can be placed in the approximate coordinates of the vector V.

The three constraint types are described in detail in sections 2.5 and 2.7.6, infra.

Variant

Until now only one certainty factor common to the x, y and z components of each constraint has been introduced but up to three different certainty factors can be assigned to these three components.

1.3.4 Encoding a surface S

In order to optimize the programming of the variant of the DSI method presented in section 2, infra, a surface composed of triangular and/or polygonal facets is encoded as a set of N atoms whose kernels (nodes) are the vertices of the triangles and/or polygons. The data and codes for each of these atoms {A(1), . . . ,A(N)} may be stored:

either in an array (see FIG. 10), or in the Atom field of memory areas linked by their Next field (see FIG. 11 ). In this case, the Next field associated with the atom number k must contain the address or the access code of the next atom numbered (k+1); the Next field associated to the last atom A(N) must contain a code indicating that the last atom has been reached.

The second solution is better if it must be possible to add and/or remove atoms easily but requires more memory than the first solution.

1.3.5 Variants

To simplify the description, the memory areas used to encode the data have been split up into fields and each of these fields has been given:

a name, a memory location.

One of ordinary skill in the art would understand that the encoding process described in this patent does not depend on the names or the order of these fields inside the memory areas used. Also, it is always possible to add new fields for particular applications without modifying the storage process.

1.4 Fitting a surface S

1.4.1 Introduction

Let S be a surface encoded according to the process described in section 1.3. There will be described hereafter a process based on this encoding to fit optimally the coordinates of some nodes {.phi..sub.1, . . . ,.phi.N} in order to:

make S as smooth as possible in relation to the global roughness criterion of the DSI method (see references [7], [8], the appendix and section 1.1.4), make S respect as much as possible the data and constraints, however accurately they may be known, concerning the shape of this surface (see sections 1.1.1 and 1.1.3).

For that purpose, we use a new variant of the DSI method described in the appendix and designed to use the encoding process described in section 1.3. Taking into account the measured or observed data, this method fits the three coordinates (.phi..sub..alpha..sup.x,.phi..sub..alpha..sup.y,.phi..sub..alpha..sup.z) of each free node .phi..sub..alpha. successively and using an iterative method and the following formula:

new value of .phi..sub..alpha..sup.x =.function..sub..alpha..sup.x (old values of .phi..sub.1.sup.x, . . . ,.phi..sub.N.sup.x), new value of .phi..sub..alpha..sup.y =.function..sub..alpha..sup.y (old values of .phi..sub.1.sup.y, . . . ,.phi..sub.N.sup.y), new value of .phi..sub..alpha..sup.z =.function..sub..alpha..sup.z (old values of .phi..sub.1.sup.z, . . . ,.phi..sub.N.sup.z).

Referring to the appendix describing the new version of the DSI method, these updating functions .function..sub..alpha..sup.x (), .function..sub..alpha..sup.y () and .function..sub..alpha..sup.z () are derived from the local form of the DSI equations at point .alpha. and have the following form (see sections 2.4.1, 2.4.2 and 2.7.3, infra): ##STR1##

The nature and the role of the weighting coefficients {v.sup..alpha. (k)} and .mu.(k) have been discussed in section 1.1.4 and are explained in article [8] and in section 2, infra. The values of the terms {.GAMMA..sub.i.sup.x.alpha.,.GAMMA..sub.i.sup.y.alpha.,.GAMMA..sub.i.sup.z.alpha. }, {.gamma..sub.i.sup.x.alpha.,.gamma..sub.i.sup.y.alpha.,.gamma..sub.i.sup.z.alpha. } and {Q.sub.i.sup.x.alpha.,Q.sub.i.sup.y.alpha.,Q.sub.i.sup.z.alpha. } depend on the types and the values of the constraints attached to the node .phi..sub..alpha. and their exact formulation is given in section 2, infra.

When the weighting coefficients {v.sup..alpha. (k)} correspond to the harmonic weighting described in section 1.1.4 and in section 2, infra, the weighting functions .function..sub..alpha..sup.x (), .function..sub..alpha..sup.y () and .function..sub..alpha..sup.z () may be simplified and take the following form: ##STR2##

Hereafter there is described a method of calculating these functions based on the encoding described in section 1.3.

1.4.2 Notation

In order to simplify the description of the computation process the following notation and definitions are used:

A block of operations is any list of operations beginning with an open curly brace "{" and ending with a closing curly brace "}". The operations in a block may be individual operations or blocks of operations.

The term working memory or working variable refers to any memory area used for storing intermediary results during a computation.

Loading the expression e in a memory m is the operation of evaluating the value of e and storing the result in the memory m. Such an operation is written as follows:

m.rarw.e

For any memory m used to store a numerical value, incrementing m by the increment i is the operation of computing the value m+i and storing the result in the memory m. Such an operation is written as follows:

m.rarw.m+i

In the above operation, the increment i may be:

a constant, or the content of a memory, or an arithmetical expression.

In the above operation, the increment i may be:

a constant, or the content of a memory, or an arithmetical expression. 1.4.3 Fitting process with harmonic weightings

When the coefficients {v.sup..alpha. (k)} are the harmonic weightings described in section 1.1.4 and in the appendix, the values .function..sub..alpha..sup.x, .function..sub..alpha..sup.y and .function..sub..alpha..sup.z taken by the functions .function..sub..alpha..sup.x (), .function..sub..alpha..sup.y () and .function..sub..alpha..sup.z () relating to the node .phi..sub..alpha. may be computed efficiently by the following method based on the encoding of S.

0) Allocate the following working memories: ##EQU7##

1) Perform the following initializations (the order does not matter): ##EQU8##

2) Let n.sub..alpha. be the number of satellites of A(.alpha.). This number can be obtained from the nb.sub.-- sat field of A(.alpha.).

3) Let .LAMBDA.(.alpha.) be the set of the satellites of the atom A(.alpha.) accessible from the Sat field of the atom A(.alpha.). For each satellite (k.epsilon..LAMBDA.(.alpha.), repeat the operations of the following block

{

3.1) Determine the number n.sub.k of satellites of A(k). This number can be obtained from the nb.sub.-- sat field of A(k).

3.2) Determine the coordinates (.phi..sub.k.sup.x,.phi..sub.k.sup.y,.phi..sub.k.sup.z) of the node .phi..sub.k of the atom A(k). These coordinates can be obtained from the Node field of A(k).

3.3) From the data in A(k) calculate or determine the value .mu.(k) of the weighting coefficient (see section 1.1.4 and section 2, infra).

3.4) Perform the following operations (in any order): ##EQU9##

3.5) Let .LAMBDA.(k) be the set of satellites of the atom A(k) accessible from the Sat field of the atom A(k). For each satellite (.beta..epsilon..LAMBDA.(k), .beta..noteq..alpha.) perform the operations of the following block:

{

3.5.2) Determine the current coordinates (.phi..sub..beta..sup.x,.phi..sub..beta..sup.y,.phi..sub..beta..sup.z) of the kernel .phi..sub..beta. of the atom A(.beta.). These coordinates can be obtained from the Node field of A(.beta.).

3.5.3) Perform the following three operations (in any order): ##EQU10## {

3.6) Perform the following three operations (in any order): ##EQU11## }

4) Perform the following operations: ##EQU12##

5) Let C(.alpha.) be the set of constraints attached to the atom A(.alpha.) that can be accessed by the Const field of this atom. For each constraint (c.epsilon.C(.alpha.)) repeat the operations of the following block:

{

5.1) Determine the type of the constraint c.

5.2) Depending on the type of the constraint c, perform the following two operations:

5.2.1) Extract the data specific to the constraint c.

5.2.2) Depending on the data selected in (5.2.1), compute the values of the following expressions (see sections 2.5 and 2.7, infra): ##EQU13##

5.3) Perform the following operations (in any order): ##EQU14## }

6) Perform the following operations (in any order): ##EQU15##

At the end of this process, the working memories .function..sub..alpha..sup.x, .function..sub..alpha..sup.y and .function..sub..alpha..sup.z contain the values of the functions .function..sub..alpha..sup.x (), .function..sub..alpha..sup.y () and .function..sub..alpha..sup.z () used to fit the values of the coordinates of the node .phi..sub..alpha. : ##EQU16##

This fitting process must be applied only to the "free" nodes, that is to say nodes that do not correspond to control nodes (see section 1.1.3).

1.4.4 Variants

Variant 1

The implementation of the fitting method based on the DSI method described in the previous section in the case of harmonic weighting coefficients may be very easily adapted to the case of any other kind of weighting. The general structure of the process remains unchanged, the only changes concerning the incrementing of the working memories allocated at step (0). How these working memories are incremented depends on the weighting adopted.

For example, if the following weighting coefficients {v.sup..alpha. (k)} are chosen: ##EQU17## then it is easy to verify that the associated fitting functions are: ##STR3##

The process used to compute the values .function..sub..alpha..sup.x, .function..sub..alpha..sup.y and .function..sub..alpha..sup.z of these functions at node .phi..sub..alpha. is then almost identical to the process described for harmonic weighting:

0) Allocate the following working memories: ##EQU18##

1) Perform the following initializations (the order does not matter): ##EQU19##

2) Determine the number n.sub..alpha. of satellites of A(.alpha.). This number can be obtained from the nb.sub.-- sat field of A(.alpha.)

3) Let .LAMBDA.(.alpha.) be the set of satellites of the atom A(.alpha.) accessible from the Sat field of the atom A(.alpha.). For each satellite (k.epsilon..LAMBDA.(.alpha.)) repeat the operations of the following block:

{

3.1) Determine the number n.sub.k of satellites of A(k). This number can be obtained from the nb.sub.-- sat field of A(k). Then perform the following operation:

n.sub.k.sup.2 .rarw.n.sub.k .times.n.sub.k

3.2) Determine the coordinates (.phi..sub.k.sup.x,.phi..sub.k.sup.y,.phi..sub.k.sup.z) of the Node .phi..sub.k of the atom A(k). These coordinates can be obtained from the Node field of A(k).

3.3) From the data in A(k) calculate or determine the value .mu.(k) of the weighting coefficient (see section 1.1.4 and section 2, infra).

3.4) Perform the following operations (in any order): ##EQU20##

3.5) Let .LAMBDA.(k) be the set of satellites of the atom A(k) accessible from the Sat field of the atom A(k). For each satellite .beta..epsilon..LAMBDA.(k), .beta..noteq..alpha. perform the operations of the following block:

{

3.5.2) Determine the current coordinates (.phi..sub..beta..sup.x,.phi..sub..beta..sup.y,.phi..sub..beta..sup.z) of the kernel .phi..sub..beta. of the atom A(.beta.). These coordinates can be obtained from the Node field of A(.beta.).

3.5.3) Perform the following three operations (in any order): ##EQU21##

}

3.6) Perform the following three operations (in any order): ##EQU22## }

4) Perform the following operations: ##EQU23##

5) Lot C(.alpha.) be the set of constraints attached to the atom A(.alpha.) accessed by the Const field of this atom. For each constraint c.epsilon.C(.alpha.) repeat the operations of the following block:

{

5.1) Determine the type of the constraint c.

5.2) Depending on the type of the constraint c, perform the following two operations:

5.2.1) Extract the data specific to the constraint c.

5.2.2) Depending on the data selected in (5.2.1), compute the values of the following expressions (see sections 2.5 and 2.7, infra): ##EQU24##

5.3) Perform the following operations (in any order): ##EQU25## }

6) Perform the following operations (in any order): ##EQU26##

The numerical parameters stored in the Info field of the atoms may be interpolated by the DSI method (see appendix: sections 2.1.2 and 2.2, infra). In this case, the process is almost identical to the process described above used to fit a surface; the only modification consists in using only one of the three functions normally dedicated to fitting the three coordinates of the nodes and to calculate it according to the interpolated parameter.

The present invention finds applications in the following fields in particular:

in geology and in geophysics, it can be used to create a surface model, for example to construct a solid three-dimensional model, on which a geophysical simulation or an underground resources exploitation simulation can be carried out; it can also be used to obtain a set of data for initializing a digital or analogue simulator of geophysical or resource exploitation phenomena; it can further be used, in conjuction with a graphics workstation, to display on a screen or to print out on paper images used to monitor simulations as above; in biology or in medicine, the invention can be used to construct three-dimensional models of organs and/or prostheses or to simulate the effects of plastic surgery; similarly, displayed or printed two-dimensional images can be used to monitor the above operations.

However, it is obvious that the invention applies more generally to modeling surfaces of any kind and in particular any natural or man-made surface. It can be of particular benefit in computer-aided design (CAD).

Finally, the present invention is in no way limited to the embodiment described above and shown in the drawings and those skilled in the art can vary and modify the embodiment as described above without departing from the scope of the invention.

NEW VERSION OF THE DSI METHOD FOR INTERACTIVE MODELING OF COMPLEX SURFACES

In Geology and Biology, a common problem is modeling complex surfaces such as interfaces between areas of different kinds or having different properties. Classical modeling techniques based on Bezier interpolations and spline functions [9] are not well suited to processing this type of heterogeneous data. A different approach is based on the DSI ("Discrete Smooth Interpolation") method [8]. In this approach, surfaces are modeled using irregular triangular facets whose vertices must be determined to allow for the largest possible quantities of heterogeneous data.

FIG. 12 shows an apparatus according to the present invention, which can be utilized to implement the different surface modelling methods of the present invention. The surface modelling apparatus, as exemplified in FIG. 12, comprises a measuring device which measures to obtain a set of geometrical data, a meshing device which meshes a surface to be modelled, a memory to store data related to and for each mode of the mesh, a calculating device which processes fitting of the surface according to the stored data, and a surface representation generator for creating a representation of the surface from fitted coordinates of each node.

2.1 Introduction: Statement of the problem

2.1.1 Surface S and related definitions

Referring to FIG. 1, let S be a surface composed of contiguous flat triangular facets embedded in the 3D Euclidian space (O.vertline.x,y,z). This surface need not be connex or closed. The following notation is used: ##EQU27##

G is a graph whose nodes are the points of .OMEGA. and this set .OMEGA. is identified with the set of integer indices {1, . . . ,N} used to number the nodes. The neighborhood N(k) of a node k.epsilon..OMEGA. iS defined as: ##EQU28##

In the following, it is assumed that these neighborhoods N(k) have the following symmetry property:

.alpha..epsilon.N(.beta.) .beta..epsilon.N(.alpha.)

The neighborhoods N(k) so defined thus generate a topology on G that can be considered as an approximate topology for S itself.

2.1.2 Problem 1: Estimating a function .phi. defined on .OMEGA.

Assume that the position of each node k.epsilon..OMEGA. iS known and let .phi.(k) be a function defined for all nodes k .epsilon..OMEGA. and known only on a subset of .OMEGA.: ##EQU29##

Generally, L is different from .OMEGA. and there is an infinity of functions .phi.(k) defined on .OMEGA. and interpolating the values {.phi.(I)=.phi..sub.I :I.epsilon.L}. The goal is to select among all such functions the one which minimizes a given criterion R*(.phi.) measuring:

the global roughness of the function .phi.(.), the discrepancy between the function .phi.(.) and some constraints and data that it should satisfy. 2.1.3 Problem 2: Fitting the surface S

In this case, only a subset of L of nodes k.epsilon..OMEGA. have a known position .phi.(k)=.phi..sub.k. ##EQU30##

As for the first problem, L is generally different from .OMEGA. and there is an infinity of surfaces S interpolating the points {.phi.(I)=.phi..sub.I :I.epsilon.L}. The goal is to select among all such surfaces the one which minimizes a given criterion R*(.phi.) measuring:

the global roughness of the surface S, the discrepancy between the surface S and some constraints and data that it must satisfy.

Allowing for the independence of the components (.phi..sup.x (k),.phi..sup.y (k),.phi..sup.z (k)) of each vector .phi.(k), it is easy to use the solution of the first problem to solve the second by simply writing:

R*(.phi.)=R*(.phi..sup.x)+R(.phi..sup.y)+R*(.phi..sup.z)

In the very particular case where the set .OMEGA. corresponds to the nodes of a regular rectangular grid laid on the surface, consideration might be given to using a Bezier, spline or Briggs type interpolation method (See [9], [1] and [3]); unfortunately these methods are not appropriate for irregular grids and cannot allow for the constraints described in sections 2.2.3, 2.5 and 2.7.2.

2.2 Interpolating a function .phi.(k) defined on .OMEGA.

This section briefly introduces the Discrete Smooth Interpolation method described in [8].

2.2.1 Notation

In the following, .phi. denotes a column matrix of size N as: ##EQU31##

The problem does not depend on the method of numbering the nodes of the grid. In order to simplify the notation, it will therefore be assumed that a permutation of the elements of the matrix .phi. has been performed in such a way that .phi. can be split into two submatrices .phi..sub.I and .phi..sub.L such that: ##EQU32##

Moreover, .parallel. . .parallel.D denotes a seminorm associated to a square (N.times.N) positive semidefinite matrix [D] in such a way that, for any column matrix X of size N:

.parallel.X.parallel..sub.D.sup.2 =X.sup.t .multidot.[D].multidot.X 2.2.2 Defining a global roughness criterion R(.phi.)

Let R(.phi..vertline.k) be a local roughness criterion defined by the following formula where {v.sup..alpha. (k)} are given positive, negative or null weighting coefficients: ##EQU33##

R(.phi..vertline.k) can be used to derive a global roughness criterion R(.phi.) defined as follows, where .mu.(k) is a given weighting function defined on .OMEGA.: ##EQU34##

R(.phi.) is completely defined by the coefficients {v.sup..alpha. (k)} and .mu.(k) and it is easy to verify (See [8]) that there is always an (N.times.N) symmetric positive semidefinite matrix [W] such that:

R(.phi..vertline.k)=.phi..sup.t .multidot.[W(k)].multidot..phi. 2.2.3 Allowing for linear constraints in a least squares sense

Consider a square (N.times.N) matrix [A.sub.i ] and an N column matrix B.sub.i defining the linear constraint:

[A.sub.i ].multidot..phi..perspectiveto.B.sub.i

For an N.times.N positive semidefinite matrix [D.sub.i ] ".perspectiveto." means:

.parallel.[A.sub.i ].multidot..phi.-B.sub.i .parallel..sub.D.sbsb.i.sup.2 as small as possible

If there are several conditions of this type, the degree of violation of these constraints can be measured by the criterion .rho.(.phi.) with: ##EQU35## 2.2.4 Solution of the problem

Among all the functions .phi.(k) defined on .OMEGA. and interpolating the data {.phi.(I)=.phi..sub.I :I.epsilon.L }, one is selected which minimizes the following criterion:

R*(.phi.)=R(.phi.)+.rho.(.phi.)

Developing R*(.phi.): ##EQU36##

The partition of the matrix .phi. induces a similar partition of the matrices [W*] and Q used to define R*(.phi.): ##EQU37##

The condition .differential.R*(.phi.)/.differential..sub.100 .sbsb.I =[0] yields the following "DSI equation" characterizing all the functions {.phi.(k):k.epsilon..OMEGA.} constituting a solution of the problem: ##EQU38##

2.3 Uniqueness of the solution of the DSI equation

Definition

The set L of nodes where .phi. is known is said to be consistent relative to R*(.phi.) if each connex component of the graph G contains at least one node belonging to L.

Theorem

If the global roughness criterion R(.phi.) is such that: ##EQU39## the DSI equation based on R(.phi.) has a unique solution.

Proof

The theorem explained above is of theoretical interest only and does not affect in any way the method described in this patent. It therefore requires no proof here.

2.4 Local form of the DSI equation

In the following, instead of solving directly the DSI matrix equation, an iterative approach avoids the computation and storage of [W*]. This iterative approach is that used in this patent.

In order to simplify the notation, it will be assumed in the following that the positive semi-definite matrices [Di] of the DSI equation are diagonal. Moreover, when linear constraints of the type [A.sub.i ].multidot..perspectiveto.B.sub.i exist, the following notation is used: ##EQU40## 2.4.1 Computing .differential.R*(.phi.)/.differential..sub..phi..sbsb..alpha.

The definition of R(.phi..vertline.k) implies: ##EQU41##

From this it can be deduced that: ##EQU42##

Allowing for the diagonal structure of [Di], it is easy to verify that the .alpha..sup.th element of .differential..parallel.[A.sub.i ].multidot..phi.-B.sub.i .parallel..sub.D.sbsb.i.sup.2 /.differential..sub..phi..sbsb..alpha. is such that: ##EQU43##

Given the definition of R*(.phi.), it is possible to deduce that: ##EQU44##

The solution .phi. is such that .differential.R*(.phi.)/.differential..sub..phi..sbsb..alpha. =0, hence the .alpha..sup.th component .phi..sub..alpha. of .phi. must satisfy the following equation: ##EQU45##

In this disclosure this equation is called the local form of the DSI equation at node .alpha..

2.4.2 Proposition for an iterative algorithm

The above local form of the DSI equation suggests a straightforward algorithm to estimate the solution .phi.. For example, at iteration (n+1) the .alpha..sup.th component .phi..sub..alpha..sup.(n+1) of the solution .phi..sup.(n+1) must satisfy the DSI(.alpha.) equation so that an iterative algorithm can be:

let I be the set of indexes of nodes where .phi..sub..alpha. is unknown, let .phi. be an initial approximate solution, while (more iterations are needed). ##EQU46##

Note that this very simple algorithm does not use explicitly the matrix [W*] occurring in the DSI equation but has to compute repeatedly products such as {v.sup..alpha. (k).multidot.v.sup..beta. (k)} which are the products used to derive [W].

In fact, if the initial approximate solution is close to the exact solution, few iterations are needed and the computation overhead becomes negligible. For example, this occurs in an interactive context where the initial solution .phi. is taken equal to the solution before some local modifications are made by the user. In this case, despite the slight increase in overhead, the local form of the DSI equation is preferable as it is much easier to use than the global form.

In order to show that the proposed algorithm actually converges, note that R*(.phi.) can be expressed as a function of .phi..sub..alpha. :

R*(.phi.)=A.multidot..phi..sub..alpha..sup.2 +B.multidot..phi..sub..alpha. +C

The coefficients A, B and C are independent of .phi..sub..alpha. and: ##EQU47##

The minimum of the global criterion R*(.phi.) is achieved for the value .phi..sub..alpha. =-B/(2A) which is precisely the value given by the DSI(.alpha.) equation. This shows that the algorithm converges because, at each step of the iterative process, the value of the positive or null function R*(.phi.) decreases.

2.5 Allowing for fuzzy data

This section introduces two types of linear constraints that will be of special importance for Geometric Modeling, as will be shown in section 2.7.6. The purpose of this section is to account for fuzzy data concerning the unknown values {.phi..sub..lambda. :.lambda..epsilon.I}. Such data is assumed to be associated with certainty factors .omega..sup.2 .epsilon.R.sup.+ which are actually positive weighting coefficients.

2.5.1 Case of isolated fuzzy data

Let .lambda..epsilon.I be a node index for which the unknown value .phi..lambda. is assumed to be close to an uncertain datum .phi..lambda. with a certainty factor equal to .omega..sub..lambda..sup.2.

.phi..sub..lambda. .perspectiveto..phi..sub..lambda. where .perspectiveto. represents the following condition:.phi..lambda.: .omega..sub..lambda..sup.2 .multidot..vertline..phi..sub..lambda..sup.* -.phi..sub..lambda. .vertline..sup.2 minimum

Such a condition is easily accounted for in the DSI method as the i.sup.th constraint with: ##EQU48##

The corresponding coefficients .GAMMA..sub.i.sup..alpha., Q.sub.i.sup..alpha. and .gamma..sub.i.sup..alpha. used in the local DSI(.alpha.) equation are thus defined by: ##EQU49## 2.5.2 Case of differential fuzzy data

Let .lambda..epsilon.I be a node index for which the unknown value .phi..lambda. is assumed to be linked to another value .phi..mu. corresponding to another index .mu..noteq..lambda.. Assume that this link is of the following form, where .DELTA..sub.80 .mu. is a given value:

(.phi..sub..mu. -.phi..sub..lambda.).perspectiveto..DELTA..sub..lambda..mu.

This relation is assumed to be fuzzy and is modeled through the following condition involving a certainty factor .omega..sub..lambda..mu..sup.2 :

.omega..sub..lambda..mu..sup.2 .multidot..vertline.(.phi..sub..mu..sup.* -.phi..sub..lambda..sup.*)-.DELTA..sub..lambda..mu. .vertline..sup.2 minimum

Such a condition can easily be accounted for in the DSI method as the i.sup.th constraint with: ##EQU50##

The corresponding coefficients .GAMMA..sub.i.sup..alpha., Q.sub.i.sup..alpha. and .gamma..sub.i.sup..alpha. used in the DSI(.alpha.) equation are thus defined by: ##EQU51## 2.5.3 Choosing the certainty factors

Consider the term M.sub..alpha..sup.* of the DSI(.alpha.) equation: ##EQU52##

In the two examples discussed above, the terms .gamma..sub.1.sup..alpha. are either equal to zero or equal to the certainty factor:

.gamma..sub.i.sup..alpha. =.omega..sub.i.sup.2

This suggests choosing for .omega..sub.i.sup.2 a given percentage P.sub.i of M.sub..alpha. :

.omega..sub.i.sup.2 =p.sub.i .multidot.M.sub..alpha. p.sub.i >0

The term M.sub..alpha..sup.* of the DSI(.alpha.) equation is then:

M.sub..alpha..sup.* =M.sub..alpha. .multidot.(1+p.sub.1 (.alpha.)+p.sub.2 (.alpha.)+ . . . ) where p.sub.i (.alpha.) is either equal to a given percentage P.sub.i or equal to zero.

2.6 Choosing the weighting coefficients

The choice of the weighting coefficients {v.sup..alpha. (k)} and {.mu.(k)} is completely free except that the {.mu.(k)} coefficients have to be positive or equal to zero. In the following an example is given of how to choose these two families of coefficients.

2.6.1 Choosing the {v.sup..alpha. (k)} harmonic weighting coefficients

Let .LAMBDA.(k) be the "orbit" of k defined as follows:

.LAMBDA.(k)=N(k)-{k}

Let .vertline..LAMBDA.(k).vertline. be the number of elements of (k). The weighting is called harmonic weighting if the coefficients {v.sup..alpha. (k)} are chosen according to the following definition: ##EQU53##

Harmonic functions have the characteristic property that they are equal at any point to their own mean on a circle centered on this point. As one can see, R(.phi..vertline.k)=0 if .phi..sub.k iS equal to the mean of the values .phi..sub..alpha. surrounding the node k, this is why, by analogy with harmonic functions, the term harmonic has been adopted for the weighting coefficients {v.sup..alpha. (k)} described in this section.

If the DSI(.alpha.) equation is developed with these coefficients, the following equation is obtained, which can be easily translated into programming language: ##EQU54## 2.6.2 Choosing the {.mu.(k)} coefficients

The coefficients {.mu.(k)} have been introduced in order to modulate locally the smoothness of the interpolation. If there is no special reason to do differently, a uniform weighting may be used:

.mu.(k)=1 k.epsilon..OMEGA.

Among all non-uniform weighting schemes, one seems to be particularly interesting if smoothness is needed in the neighborhood of the data: ##EQU55##

2.7 Application to Geometric Modeling

This section shows how to use the DSI method to estimate the triangulated surface S itself for which it will be assumed that partial data is available composed of:

the exact location of some apices of the triangles, the approximate location of some apices of the triangles, vectorial constraints specifying the shape of S, etc.

To this end, .phi..sub..alpha. is defined as the current vector joining the origin of R.sup.3 to the current node of .OMEGA. and .phi..sup.x, .phi..sup.y and .phi..sup.z denote its components on the {ox,oy,oz} orthonormal basis of R.sup.3 :

.phi.=(.phi..sup.x,.phi..sup.y,.phi..sup.z) 2.7.1 Defining a global roughness criterion R(.phi.)

As usual, the set of nodes .OMEGA. is split into two subsets I and L such that: ##EQU56##

The nodes I.epsilon.L whose location .phi..sub.I .epsilon.R.sup.3 is known will be called control nodes and the aim is to determine the location of the remaining points i.epsilon.I in such a way that the following local criterion R(.phi..vertline.k) based on a given set of weighting coefficients {v.sup..alpha. (k)} is as small as possible: ##EQU57##

This criterion constitutes the vectorial form of the local DSI criterion. The associated global vectorial DSI criterion R(.phi.) is defined as follows where {.mu.(k)} are given non-negative weighting coefficients:

R(.phi.)=.SIGMA..mu.(k).multidot.R(.phi..vertline.k)

Defining R(.phi..sup.x), R(.phi..sup.y) and R(.phi..sup.z) as previously: ##EQU58##

Pythagoras' theorem yields:

R(.phi.)=R(.phi..sup.x)+R(.phi..sup.y)+R(.phi..sup.z) 2.7.2 Allowing for linear constraints in a least squares sense

Consider the following matrices sized so that N is the total number of nodes of the set .OMEGA.: ##EQU59##

To have the components .phi..sup.x,.phi..sup.y and .phi..sup.z satisfy the following constraints

[A.sub.i.sbsb.x.sup.x ].multidot..phi..sup.x .perspectiveto.B.sub.i.sbsb.x.sup.x [A.sub.i.sbsb.y.sup.y ].multidot..phi..sup.y .perspectiveto.B.sub.i.sbsb.y.sup.y [A.sub.i.sbsb.z.sup.z ].multidot..phi..sup.z .perspectiveto.B.sub.i.sbsb.z.sup.z the following criterion is introduced: ##EQU60## 2.7.3 Solution of the problem

Among all the surfaces S interpolating the data {.phi.(I)=.phi..sub.I :I.epsilon.L }, the aim is to select the one minimizing the following criterion:

R*(.phi.)=R(.phi.)+.rho.(.phi.)

It is easy to verify that this equation can be rewritten as follows: ##EQU61##

Allowing for independence of the components .phi..sup.x,.phi..sup.y and .phi..sup.z : ##EQU62##

Minimizing the vectorial DSI criterion is thus equivalent to minimizing the corresponding three DSI criteria applied to the three components of the current vector .phi..

2.7.4 Note

If the surface S is composed of several disjoint patches the theorem explained in section 2.3 ensures that there exists a unique solution to the problem if at least one triangle apex within each patch has a known position.

2.7.5 Interactive use

In Geometric Modeling, the goal is to define interactively the geometric shape of the graph G associated with the set of nodes .OMEGA.. To this end, it is assumed that an initial shape is known and the user moves some nodes and changes interactively some constraints (control nodes and fuzzy data); these modifications are then propagated to the nodes i.epsilon.I through the vectorial form of the DSI method.

As explained previously, this problem is broken down into three sub-problems corresponding to the three coordinates (x,y,z) and the best way to solve these three problems is to use the iterative algorithm associated with the local form of the corresponding DSI equations.

Between steps of the modeling process, not only can the position of the control nodes be modified but also the subset L itself can be changed if necessary; for example, to apply DSI only to a subset .GAMMA. of nodes, it is necessary to choose L in such a way that .GAMMA. .GAMMA. where .GAMMA. stands for the complementary set of .GAMMA. in .OMEGA..

Moreover, because the problem is broken down into sub-problems solved independently and corresponding to the three coordinates, different subsets L.sub.x, L.sub.y and L.sub.z can be defined for each coordinate if necessary; for example, if the (x,y) coordinates must not be modified it is sufficient for L.sub.x and L.sub.y to be equal to the whole set .OMEGA..

2.7.6 Allowing for fuzzy geometrical data

Initially, the DSI method was designed to model complex surfaces encountered in the field of natural sciences, for example in geology and biology. In this case, there is generally much imprecise data to take into account such as: ##EQU63##

Fuzzy control nodes

By definition, a "fuzzy control node" is any node .lambda..epsilon.I whose position .phi..lambda. is known as being approximately equal to a given vector .phi..lambda. with a given degree of certainty .omega..sub..lambda..sup.2. Referring to section 2.5.1, such fuzzy data can be taken into account with the DSI method by introducing the following constraints in the R*(.phi..sup.x), R*(.phi..sup.y) and R*(.phi..sup.z) criteria: ##EQU64##

Fuzzy vectorial constraints

In many applications it is necessary to control the vector (.phi..sub..mu. -.phi..sub..lambda.) joining two nodes .phi..sub..lambda. and .phi..sub..mu. of a surface. By definition, such a constraint is called a "fuzzy vectorial constraint" if it must be satisfied with a given degree of certainty .omega..sub..lambda..mu..sup.2. Some examples of constraints of this type are shown in FIG. 2.

In practice, at least one of the two points .phi..sub..lambda. and .phi..sub..mu. has an unknown position, say the one corresponding to .lambda..epsilon.I, and two important cases have to be considered:

The first case corresponds to FIGS. 2a and 2b where the shape of an object must be controlled in such a way that (.phi..sub..mu. -.phi..sub..lambda.) is equal to a given vector .DELTA..sub..lambda..mu..

The second case corresponds to FIG. 2c where the shape of an object must be controlled in such a way that (.phi..sub..mu. -.phi..sub..lambda.) is orthogonal to a given vector V. This situation occurs when the vector normal to a surface is to be controlled; in this case, if V is the normal vector at node .lambda. of the surface, it is wise to control its orthogonality with all the vectors (.phi..sub..mu. -.phi..sub..lambda.) for all .mu..epsilon.N(.lambda.).

Referring to section 2.5.2, the first case is easily taken into account in the DSI method by introducing the following constraints in the R*(.phi..sup.x), R*(.phi..sup.y) and R*(.phi..sup.z) criteria with a weighting coefficient equal to .omega..sub..lambda..mu..sup.2 : ##EQU65##

For the second case, the problem is not very different because: ##EQU66##

In this case, the values (.phi..sub..lambda..sup.x,.phi..sub..lambda..sup.y,.phi..sub..lambda..sup.z) and (.phi..sub..mu..sup.x,.phi..sub..mu..sup.y,.phi..sub..mu..sup.z) on the righthand side of the equations can be set to the corresponding values inferred at the previous step of the iterative process; if

.DELTA..sub..lambda..mu..sup.x =-{V.sup.y .multidot.(.phi..sub..mu..sup.y -.phi..sub..lambda..sup.y)+V.sup.z .multidot.(.phi..sub..mu..sup.z -.phi..sub..lambda..sup.z)}/V.sup.x .DELTA..sub..lambda..mu..sup.y =-{V.sup.z .multidot.(.phi..sub..mu..sup.z -.phi..sub..lambda..sup.z)+V.sup.x .multidot.(.phi..sub..mu..sup.x -.phi..sub..lambda..sup.x)}/V.sup.y .DELTA..sub..lambda..mu..sup.z =-{V.sup.x .multidot.(.phi..sub..mu..sup.x -.phi..sub..lambda..sup.x)+V.sup.y .multidot.(.phi..sub..mu..sup.y -.phi..sub..lambda..sup.y)}/V.sup.z the second case becomes identical to the first case.

Note

In allowing for fuzzy data as described above only one certainty factor .omega..sub..lambda..mu..sup.2 is used for all three components of the constraints. In practice this is not mandatory and, because of the independence of the three components, it is always possible to use three different certainty factors .omega..sub..lambda..mu..sup.x2, .omega..sub..lambda..mu..sup.y2 and .omega..sub..lambda..mu..sup.z2 corresponding to the three components.

2.8 Generalization

Before concluding, some straightforward generalizations of the approach are worth mentioning:

The DSI method can be used on any kind of polygonal facets. In practice, triangles are recommended since they are the simplest polygons and are easy to handle in a computer program.

The DSI method can be used to estimate polygonal curves in a 3D space; in this case, triangles are replaced by segments.

2.9 Conclusion

The local form of the DSI equation provides a simple and powerful tool for modeling complex surfaces encountered in Geology and Biology, for example. In particular:

There are no restrictions as to the mesh used to model a surface and it is possible to uses automatic algorithms to build them; for example, if the mesh is composed of triangles, it is possible to use algorithms derived from the Delaunay method (See [2] and [5]). Moreover, the size of the mesh can easily be fitted locally to the complexity of the surface.

The DSI method can cope with large quantities of precise or fuzzy data. This may be particularly useful in Natural Sciences where the aim is not to generate aesthetic surfaces but to fit precise and fuzzy data consistently.

The algorithm derived from local DSI equation is fast enough to allow interactive use.

Nothing has been said about the representation of the facets. Depending on the application, plane facets may suffice. However, facets can if necessary be represented by non-planar surface patches interpolating the nodes of the mesh (See [4]).

The invention being thus described, it will be obvious that the same may be varied in many ways. Such variations are not to be regarded as a departure from the spirit and scope of the invention, and all such modifications as would be obvious to one skilled in the art are intended to be included within the scope of the following claims.

Claims

1. A process for modeling a surface representing an interface between two areas with different properties in a three-dimensional body, comprising the steps of:

measuring to obtain a set of geometrical data, from said interface of said three-dimensional body, relating to the surface and associated with respective points on said surface;

meshing the surface so that all said points are a subset of nodes of a mesh;

storing at a specific memory address for each node of the mesh, the following data:

first, second, and third coordinates,

a number of satellite nodes,

satellite node address data providing access to specific coordinates of said satellite nodes and consequently to satellite node related data relating thereto, and

geometrical data associated with said each node of the mesh;

for each node of the mesh, defining a local roughness index derived from a weighted sum of the coordinates of the node and of its satellites;

defining a global roughness index obtained by summing the local roughness indices associated with the nodes, and a global violation index of said geometrical data, to thereby generate a sum of said global roughness index and said global violation index;

fitting fitted coordinates of each node processed by an iterative method in which for each step of an iteration there is added a weighted combination of the coordinates of the node and of the satellites of said node and a combination of the geometrical data associated with said node, to thereby minimize said sum of the global roughness index and the global violation index; and

creating a representation of the surface from the fitted coordinates of each node.

2. A process according to claim 1, wherein:

the surface is meshed using triangles.

3. A process according to claim 1 wherein:

the fitting step comprises three separate sub-steps respectively using the first, second, and third coordinates of the nodes and the coordinates of the satellite nodes corresponding thereto, respectively.

4. A process according to claim 1, wherein:

the geometrical data includes vector data between two given nodes.

5. A process according to claim 1, wherein:

the geometrical data includes data of a vector normal to the surface to be modeled.

6. A process according to any one of claims 4 to 5, wherein:

at least one geometrical datum is associated with a coefficient representing a certainty factor by which said at least one datum is known, the global violation index is a function of said certainty factors.

7. A process according to claim 1, further comprising the steps of measuring properties of the three-dimensional body in a region of at least one node of the mesh and storing the measured properties at addresses specific to the at least one node so as to obtain complementary data of the three-dimensional body.

8. A process according to claim 1, wherein:

for each local roughness index, a ratio between a first and a second weight of a corresponding node equals a negative of the number of satellite nodes.

9. A process according to claim 8, wherein:

said weighted combination uses substantially the same weighting as the ratio.

10. A process according to claim 1, wherein:

the summing of the local roughness indices is a weighted summing.

11. A process according to claim 10, wherein:

said weighted combination uses the weighting associated with the local roughness indices.

12. A process according to claim 1, wherein:

the three-dimensional body is a geological formation.

13. A process according to claim 1, wherein:

the surface represents variations of a property of a two-dimensional body in a vicinity of a plane intersecting said body, said property being represented by a remaining dimension, and a representation of the surface is used to optimize visualization of underground resources of said three-dimensional body, said three dimensional body being a geological formation.

14. A process according to claim 1, wherein:

the representation of the surface is a graphical representation on a flat medium.

15. A method according to claim 1, wherein:

the representation of the surface is a three-dimensional model.

16. A device for modeling a surface representing an interface between two areas with different properties in a three-dimensional body, comprising:

means for measuring to obtain a set of geometrical data, from said interface of said three-dimensional body, relating to the surface and associated with respective points on said surface;

means for meshing the surface so that all said points are a subset of nodes of a mesh;

means for storing at a specific memory address for each node of the mesh, the following data:

first, second, and third coordinates,

a number of satellite nodes,

satellite node address data providing access to specific coordinates of said satellite nodes and consequently to satellite node related data relating thereto, and

geometrical data associated with said each node of the mesh;

calculation means for fitting fitted coordinates of each node processed by an iterative method in which for each step of an iteration there is added a weighted combination of the coordinates of the node and of the satellites of said node and a combination of the geometrical data associated with said node, to thereby minimize a sum of a global roughness index obtained by summing local roughness indices associated with the nodes, each local roughness index derived from a weighted sum of the coordinates of the node and the coordinates of its satellites, and of a global violation index of said geometrical data; and

means for creating a representation of the surface from the fitted coordinates of each node.