Camera System and Methods for Extracting 3D Model of Viewed Object

This invention relates to vision systems and, more particularly, to a novel system and methods for extracting a 3D model of an object present in a plurality of views, from at least one camera, and describable by geometric polyhedrons subject to constrained nonlinear optimization.

It has been the objective of considerable work in the art to develop and implement a complete system with which a human user can reconstruct the 3-dimensional geometry of a scene from images. Emphasis has been placed on the development of a practical, robust interactive system in which the operator guides the object selection process and the system computer provides an optimal numerical result. Although there are many factors (e.g. surface reflectance, texture, lighting and the like) that determine how a scene will appear in an image, the most directly exploitable factor is the projective geometric relationship between scene points and image points. It is highly desirable to allow the human operator to perform the high level perceptual functions of surmising the general layout of the scene and designating the locations of scene features, functions for which a human is extremely adept, and have the computer provide the optimization to recover the scene geometry of maximum likelihood. Measurements made by the human operator are considered to be noisy observations, as image resolution will limit the accuracy to which image measurements can be made. The projective relationship is described by the equations of the camera model. The recovery of camera model parameters (camera positions, orientations, and focal lengths) will necessarily be part of the process of reconstructing the scene, as this information is normally unavailable or unreliable a priori.

Besides image measurements, the other primary means of reconstructing scene geometry will be through the specification of model geometric constraints. Common constraints arising in man-made structures are perpendicularity and parallelism. Constraints due to the polygonal mesh scheme used to model surfaces must also be taken into consideration in the solution process. It is also desirable to include a provision for constraints, to provide a means for solving what otherwise would be an underspecified problem, due to insufficient information visible in the images. Also, constraints are a means for allowing a human operator to use the system by exercising control over the models generated by the system. The system must be established (`programmed`) to guard against constraint-conflict problems, especially with a human operator in the loop.

Many other forms of prior scene reconstruction methods have been developed within the field of photogrammetry, which is defined, in The Manual of Photogrammetry, as "the science or art of obtaining reliable measurements by means of photographs." It is known that one of the fundamental problems of photogrammetry is the camera modeling problem: for a given photograph, determine a mathematical model which describes how the image was formed, i.e. how points in the three-dimensional scene are projected onto a sensor film. In photogrammetry, the scene model is normally a set of 3-D points. The basic measurements which are provided from which to determine the camera model are: 1) the 2-D image coordinates of a set of relevant 3-D points; 2) the 3-D coordinates of some of these points, which are called control points; and 3) auxiliary data such as known camera parameters (such as focal length, pixel spacing in the image sensor, and the like). The control point positions are obtained from geodetic surveys or from topographic maps. The procedure to estimate the positions of all 3-D points as well as all the parameters of the camera model is to apply the principle of maximum likelihood. All measurements (image coordinates, control point positions, measured camera parameters) are assumed to contain a random error component, modeled with the normal distribution. The estimate for the 3-D point positions and camera parameters is that which corresponds to the most likely set of random errors. Finding this solution requires finding a least squares solution to a set of nonlinear equations. It is desirable to provide a system and method to extend the photogrammetry method to a more general and complex scene model.

A single view of a scene, being a two-dimensional projection of a three-dimensional scene, does not usually provide enough information to uniquely reconstruct the scene geometry. The situations with multiple views are much better, even if the viewpoints of the cameras are unknown. In computer vision, the reconstruction problem from multiple views has been called the structure-from-motion problem. Here, the scene is modeled as a collection of 3-D points which undergo rigid body motion. Multiple views of the scene are acquired at different times. Because of the rigid body assumption, the information collected is equivalent to having multiple simultaneous views of a stationary scene. The structure-from-motion problem very closely parallels the photogrammetry problem, with the major difference being that in structure-from-motion, nothing is assumed to be known about the scene (i.e. there are no control points). One of the major concerns in structure-from-motion has been determining the conditions under which a structure problem is solvable. The other major concern in structure-from-motion problems is the development of methods for recovering structure. But, as noted in the May 1989 special issue on visual motion in Vol. 11, No. 5 of the "IEEE Transactions on Pattern Analysis and Machine Intelligence", solution techniques have hitherto been ad hoc, i.e. not founded on sound principles of estimation theory. Also noted was the fact that the structure-from-motion problem is ill-conditioned and requires precise location of points in the images. Solutions are difficult to obtain, especially given just two perspective views of a small number of points. It is therefore desirable to remove these types of problems from the scene reconstruction techniques of the present invention, preferably by the inclusion of geometric constraints.

One version of the reconstruction problem which has received considerable attention is the extraction of solid models from orthographic views, i.e. from engineering drawings. Sutherland implemented a system where a user would indicate 3-D point position on a large digitizing table by designating the point in two perpendicular, orthographic views with a digitizing pen. The calibration of the view with respect to a single coordinate frame required designating "exemplary" points, i.e. known control points. The 3-D point data was then used to build up a polyhedral model. This system is based on a three step calibrate-measure-model procedure. Step one is to calibrate the views using known control points. Step two is to use corresponding points in multiple views to measure 3-D point positions from the stereo view-ray intersections. Step three is to use the 3-D point data to build the 3-D model of the object. This procedure surrenders accuracy because estimates made at each stage are not updated given further measurements and introduction of constraints at later stages. This loss of accuracy may be adequate when working with engineering drawings, but is perhaps inadequate for lower resolution imagery with less than perfectly defined edges. It is desirable to progress beyond this "open loop" approach, as by use of the single overall optimization techniques of the present invention.

Markowsky and Wesley have described in detail their implementation of the third stage of the calibrate-measure-model process. They use a "fleshing" algorithm to reconstruct polyhedral objects from at least two projections of the wireframe of the object. The topology of the object is automatically discovered by searching for all possible topologies that could generate the given projections. However, the algorithm was based on the assumption that the directions of the parallel projections (camera models) are known, and that the projections are extracted without measurement errors.

Another approach to scene reconstruction is an "open loop" approach based on parameterized models. One example is the SRI Cartographic Modeling Environment, which allows the user to interactively manipulate the shape and placement of a 3-D model within the overall 3-D scene model so that the projection of the model aligns with an object in the image. The 3-D models are generally simple in that they have only a few parameters available for manipulation, although with models such as the superquadric, fairly complex shapes can be obtained. The major limitation with this approach is that a human operator cannot effectively adjust all the variables to align the model with more than a single view, hence multiple views do not add precision to the models. Camera parameters computed prior to model extraction are not improved as models are generated. If models fail to align properly in multiple views, one view is selected as the "correct" view, and all others are adjusted to minimize inconsistencies.

Constraint-based solid modeling as a mechanical design tool has been a focus of research for several years at the MIT Computer Aided Design Lab. In the first system developed, a part would be designed through the two-step process of defining the topology of the part using a boundary representation and then using three-dimensional constraints to position characteristic points of the part. Given this definition of a part, changes, could be made in the dimensions of constraints, and the effect on the total part would be found using Newton's method to resolve the constraint equations. This basic technique of working with a set of constraint equations was refined by applying graph matching algorithms to match constraint equations to variables. The matching created a tree-like graph which simplified calculations because variables could be solved for sequentially, except for variables contained within loops in the graph (called "strong components") which would still require simultaneous solution.

This whole approach to part design presented several problems, probably the most serious being the difficulty of providing a sufficient set of constraints to constrain a part uniquely without overspecifying the system or generating an inconsistent set of constraints. It is desirable to avoid these problems, to the greatest degree possible, in any new system and techniques for extraction of models from visual data, because constraints alone will not determine an object's dimensions. Primarily, observations from images will determine dimensions. Constraints will not have to totally specify the objects, therefore avoiding inconsistent constraints will be simpler.

Various investigations have been conducted of automatic interpretation of line drawing images as 3-dimensional polyhedral solids. The input to these systems is a set of points on the drawing and a list of lines which connect points. The lines represent projected visible edges of a set of 3-D polyhedrons. The main contribution of the line interpretation systems is the provision of a means for establishing a correct labeling of lines, but with a major flaw in that error sources are not modeled and no advantage is taken of redundant measurements (redundant measurements are detected and then ignored).

In accordance with the invention, a system for extracting a three-dimensional model of an object viewed in a selected scene, includes a plurality of viewing means for acquiring scene information, each from a different spatial-temporal aspect of the total scene; means for digitizing the scene information, for at least temporary storage in memory means associated with computer means for data processing; and an interface means by which a human operator can view the stored scene information and (1) designate observable vertices of the scene objects, (2) create a non-numeric structure graph representation of the scene (polygonal mesh), and (3) specify geometric constraints (point, line direction, line length and/or coincidence). The information generated by the user is placed in storage. The processing means retrieves structure graph and constraint data from storage and uses that data to synthesize a dependency graph, thereby eliminating the constraints. A set of graph unknowns is initialized to create a trial solution of the legal dependency graph. Thereafter, a maximum-likelihood optimization method, such as Gauss-Newton iteration, is used to vary the set of unknowns to improve the solution until the least error is present. The polygonal-mesh data found for the least-squares-error solution is then properly formatted and made available at a system output for subsequent use as needed.

In a presently preferred embodiment, the Gauss-Newton iteration can be replaced by a Levenberg-Marquardt iteration, to provide more robust convergence when a poor initial estimation occurs at the beginning of the optimization sequence.

Accordingly, it is an object of the present invention to provide a novel system and methods for extraction of a 3D model, from stored scene information, for an object viewable in the stored scene.

This and other objects of the present invention will now become apparent to those skilled in the art, upon reading the following detailed description of a presently preferred embodiment, when considered in conjunction with the associated drawings.

FIG. 1 is a schematic block diagram of a computer vision system, in accordance with the invention, for extracting a 3D model of an object from scene information;

FIGS. 1a-1d are, respectively, a perspective view of a scene, and several 2D views of the 3D scene, illustrating a portion of the problem;

FIG. 2 is a view of a simple object, considered as an example in illustrating the present invention;

FIG. 3 is a flow graph illustrating the overall methodology of the present invention;

FIG. 4 is a 3D representation of one face of an object to be extracted by the method of the present invention, and illustrating several conventions used therewith;

FIG. 5 is a dependency graph for the building shown in FIG. 2;

FIGS. 6a-6d are examples of various 3D models extracted by use of the present invention; and

FIG. 7 is an orthographic projection model of focal and camera-local planes, and useful in understanding certain simplifying conventions used in the present invention.

Referring initially to FIG. 1, a model of an object 10 is to be acquired by a system 11, having at least one camera means 12. The particular object, in a set of a plurality of objects viewable in any one scene, is selected by a human operator (not shown) interfacing with the system via a viewing means 14. The system needs data obtained by viewing the scene from a plurality of temporal-spatial locations, i.e. either a plurality of cameras view the scene at the same time from different directions, or data is gathered from the scene over a period of time as the scenic objects move, translationally and or rotationally, past a single camera. Preferably, a plurality m of camera means 12a, . . . ,12m are used, with the video output of, and pointing-direction information from, each camera being provided to a separate associated input 16a, . . . ,16m of a source selection and digitizer means 16-1 of the system signal processing electronics means 16. Advantageously, means 16-1 includes source selection multiplexers, analog video-to-digital data conversion means and the like, as well known to the art. The digital data converted from each camera is separately stored in a mass data storage memory means 16-2, under control of a data-processing system computer means 16-3, which may include a dedicated constrained least-squares computational processor as known to the art (perhaps as in U.S. Pat. No. 5,018,065, and the like). The computer includes random-access memory (RAM) means 16-4 for temporary working storage, and read-only memory (ROM) means 16-5, in which at least a portion of the operational program for means 16-3 (or all of means 16) may be stored. Selection of the portion(s) or segment(s) of the program, or operation, can be made via an operator interface keyboard 18-1 and/or a light pen 18-2 and the like I/O devices. Advantageously, the operator uses the light pen 18-2 to select a point on display 14 which determines the particular displayed object 10' to be extracted from the scenic totality. The digital data characterizing the extracted object can be used to display a graphic representation of the extracted object on means 14, if desired, and can also be made available at a system I/O means 16-6 for coupling to an output device (display, plotter, numerically-controlled machining means and the like) which may use the data to provide an object replica and the like, end uses.

As seen in FIG. 1a, the plurality of spatial-temporal aspects can be acquired as a number of views of a scene are taken, either simultaneously from several different directions, or in time sequence as the scene objects move (translationally and/or rotationally) past the camera(s). Here, a first object 20 is resting upon a surface of a second object 22, while three different camera means observe the objects from the directions of view arrow 1b, 1c and 1d, respectively. The scene information available from the individual camera means are shown in view FIGS. 1b-1d, along with a representation of the 3D vectors of Cartesian-coordinate space for the direction data associated with that view.

In accordance with one aspect of the present invention, scene reconstruction is formulated as a maximum-likelihood-estimation problem, which can be solved with a constrained, nonlinear least-squares optimization. My scene reconstruction formulation is somewhat parallel to the photogrammetric approach, but with a considerably more general scene model in which a plurality of m different images each view parts of the scene to be modeled from different vantage points. The positions and focal lengths of the camera means 12a-12m may be unknown. The problem is to recover the scene geometry from the positions of features in the stored images, each containing a scene model holding some number of objects. The position of each object is represented by a coordinate frame consisting of a rotation matrix and translation vector. Objects themselves are modeled as a constrained polygonal mesh, consisting of vertices, lines and planar faces. The polygonal mesh represents the visible side of object surfaces. Each object consists of a number n of vertices v.sub.i, where i=1, . . . ,n (i.e. each vertex is represented as a 3-D vector). Each face of an object is defined by a sequence of vertices clockwise around the face when viewed from the visible side. Lines connect adjacent vertices on the boundary of a face. Constraints can be placed on line directions or face normals by use of direction vectors.

The information describing the relationships of all the geometric entities is stored in a scene structure graph, which is a list of objects; the objects themselves contain lists of vertices, faces, and lines; faces contain lists of vertices and lines; etc. The structure graph is the topological description of the scene.

The parameters which instantiate the topology into completely specified models are the scene parameters. Each geometric entity in the model is defined by some number of parameters, which may, for convenience, exceed the minimum number of parameters or degrees of freedom (DOF) needed to specify the entity. Each of a plurality n of vertices, as noted above, is defined by a 3D vector v.sub.i, for 1.ltoreq.i.ltoreq.n. Lines are defined by a line direction vector a and a vector offset o. Points x lie on the line if the equation x=ka+o is satisfied for some scalar k. Faces are defined by a face normal vector a and any point on the plane by a vector p. Points satisfying a.multidot.(x-p)=0 lie on the plane. Determination of the scene variables (the totality of all scene parameters together) is the crux of the reconstruction problem.

In accordance with another aspect of my invention, constraints are placed upon scene geometry. Constraints implicitly arise from the definition of an object model as a planar mesh, because all vertices in a face must lie in the plane of the face. This implicit constraint is, expressed by the vector equation a.multidot.(x-p)=0 for each vertex in a face.

Additional, explicit constraints on the scene geometry may be needed; explicit, operator-determined constraints will be necessary in those situations where the camera model and scene geometry are underspecified. For example, constraints may be used to fill in parts of a scene not viewed in any image. Explicit constraints may also be provided because the human user in the display-control-and-definition loop will often desire to force scene models to satisfy particular conditions (e.g. verticality of walls and the like). My presently preferred embodiment of method and system provides for use of any one, or a combination, of four types of constraints: 1) point constraint--the user directly specifies the coordinates of a vertex of any scene object; 2) direction constraint--line directions or face normals are constrained to lie along a particular direction, which can be either fixed or left free to vary; 3) length constraint--the length of a line is fixedly established; and 4) coincidence constraint--a vertex, line or face of a selected object is constrained to be coincident with another vertex, line, or face, generally of another object in the same scene.

Referring now to FIG. 2, an illustration of the use of explicit constraints is shown. A simple house facade model 24 consisting of four polygons 26a-26d, defined by eight vertices 27-1 through 27-8 and lines 28a-28k therebetween. A constraint file "object house" can be constructed to describe the topology; such a constraint file for this simple model, using vertices 1-8 and faces defined by clockwise-arranged vertex sets, might be:

The constraints force the wall edges to be vertical and establish a coordinate frame attached to the house, at vertex 27-1 (i.e. vertex number 1). In my current system, the user generates the constraint file via the graphical user interface 14, without the necessity for typing the file into the system.

My system methodology operates on a scene with multiple objects by causing the user to select one of the object coordinate frames as the base frame for the overall scene. The computer means then determines the position of all other object frames relative to the base frame. In landscape situations where maps of a geographic area are available, the user may wish to have the scene defined with respect to a geocentric or local observer frame defined at a nearby point in the "world". The user simply creates a "world" object with a set of vertices corresponding to known points in the world. The known vertex positions are entered as point constraints, within the object's coordinate system. By selecting the "world" object frame as the base frame (and designating the positions of these control points in the images), the solutional method will then output result data for all other objects with respect to the "world" frame. The user can work within a world frame if a map is available, or within a local frame attached to one object if nothing but the images is available, all without modifying the solutional program in means 16. Without control points, however, the user must enter at least one length constraint to provide scale.

The graphical user interface can be a workstation-based unit enabling the user to conveniently create the informational input to the scene geometry solver. This type of interface allows the user to designate, via some movable means (a mouse, a light pen and the like), the positions of at least one selected vertex in the images displayed on the workstation screen (acting as means 16). The user can then directly enter the topology information by selecting sequences of vertices to define each polygonal face of the model. Constraints are similarly input by use of pull-down menus and the user directly indicating particular features in the scene structure by selecting, with the movable selection means, designated vertices on the display screen.

In accordance with another principle of my invention, the object-extraction problem is formulated by collecting all of the parameters associated with all the geometric entities in the scene into only one k-dimensional vector x which represents the scene model. Constraints on the scene model are expressed by at least one nonlinear equation, and typically a plurality of such nonlinear equations, all which can be collected into one vector equation: h(x)=0.

The positions of the vertices can be extracted from the scene parameter vector x. Let the position of the i-th vertex be denoted by the vector v.sub.i (x). These 3D points are viewed, in the aforementioned plurality of spatial-temporal aspects, by one or more cameras. The remaining variables in the system are the unknown camera parameters c.sub.j, for 1.ltoreq.j.ltoreq.m cameras. The predicted image position of the i-th vertex in the j-th image is obtained from the camera parameters by applying a 3D to 2D camera projection function, p(v.sub.i, c.sub.j). The form of this function will depend on the type of camera used (frame, panoramic, etc.).

For some combinations of vertex and image, an observation of the vertex in the image will be available. The observed 2D position of vertex i in image j will be denoted, i.sup.o.sub.ij, and treated as a noisy observation with standard deviation .sigma..sub.ij in both the horizontal and vertical directions. The set of ij pairs for which an observation is available will be denoted by the set D. For some vertices, a direct observation of its 3D position may be available, generally from a map or blueprint. The observed 3D position of the i-th vertex will be denoted v.sup.o.sub.i having standard deviation .sigma..sub.i in each of x, y and z. The set of individual (i), for each of which a 3D control point is available, will be noted by the set C.

I obtain the vector g.sup.o by concatenating all the observations i.sup.o.sub.ij i, for ij.epsilon.D, and v.sup.o.sub.i, for i.epsilon.C. The likelihood function is the probability of observing the observations given the unknowns to be estimated. Assuming zero mean, independent, Gaussian random errors, the likelihood function p(g.sup.o .vertline.x,c.sub.1, . . . ,c.sub.m) is: ##EQU1## I have found that the scene variables and camera parameters which satisfy the constraints and maximize the likelihood function is equivalent to the following optimization: ##EQU2##

This formulation combines all observations and known constraints into a single overall optimization that obtains the scene geometry and camera models corresponding to the most likely set of noise samples.

Those skilled in the art will now appreciate that the optimization described by equations (1) and (2) can be approached in various ways. The method of my invention, as shown in the flow graph of FIG. 3, first retrieves the image data from storage and, in step 30, displays the image data for the user. The user selects at least one vertex of the desired object to be extracted from the scene and, interfacing via the workstation I presently use for means 16, interactively inputs, in step 32, the scene topology data and then inputs, in step 34, data for any observed vertices found in initial processing. In step 36, the initial data is displayed for prompting the user to input constraint information. Then all of the available data is used, in step 38, to convert the constrained optimization to an unconstrained optimization by eliminating constraints simultaneous with a reduction in the number of unknowns. Parameters are grouped into constrained geometric elements with dependencies between elements recorded in an acyclic dependency graph, configured in step 40, and more fully described hereinbelow in solution section A. This dependency graph expresses the polygonal mesh as a kind of kinematic structure, where each variable (degree of freedom, or DOF) is similar to a joint. The process also checks on the consistency of constraints, because inconsistent constraints preclude an acyclic graph. The existence of a dependency graph also guarantees that the combination of constraints and observations is sufficient to uniquely specify the scenic system.

Once a legal dependency graph has been created, the system proceeds toward a solution by first, in step 42, initializing all the unknowns in a bootstrap procedure described in greater detail hereinbelow in section B. The final solution is obtained in step 44, using a Gauss-Newton iteration to improve the dependency graph, as described in section C infra. The data for the final dependency graph is then properly formatted in step 46, in accordance with a preselected format standard, and the formatted data is output, in step 48, to user-selected output devices (e.g. the display means, a plotter, a numerically-controlled machine and the like) for subsequent usage as desired.

In my novel method, constraints are eliminated, in the problem restructuring step 46, by grouping the unknowns associated with the geometry into subsets associated with the basic geometric entities of the model: coordinate frames, vertex points, lines, face planes, and direction vectors. The constraints (point, direction, length, coincidence) each define a relationship between two elements. For example, if a face contains a vertex, then the face plane and the vertex point must intersect. To eliminate this coincidence constraint, the symmetry is broken and one of the elements is placed as a constraint on the other. In the face/vertex example, either the face is free and the vertex must lie in the plane of the face, or the vertex is free and the face must intersect the vertex. The constrained element has reduced degrees of freedom (fewer variables) due to the constraint. Continuing the example and referring to its illustration in FIG. 4, if the face is constrained by the vertex, then the face is no longer represented by a 3-DOF free plane but by a 2-DOF plane 50 which is attached to a fixed point 52 in space, with DOF vectors 54a and 54b. In effect, the constraint disappears from the overall system and is submerged into the constrained geometric element.

As discussed briefly hereinabove, eliminating all constraints requires finding an acyclic graph, called a dependency graph, which defines an evaluation ordering of all elements. The ordering is such that if one element is dependent on other elements, then that one element must appear later in the evaluation ordering than all of the elements it depends on. In addition, individual geometric elements must not be overconstrained. For example, a plane cannot be constrained by four points. The allowable types of constrained geometric elements are listed in the following Table:

When generating the dependency graph, the assignment of node types to geometric entities is restricted to insure that during the subsequent solution process, initializing unknowns will be straightforward. For example, a vertex which is not observed in any of the images of a problem set is restricted to node types which allow the position of the vertex to be deduced solely from its predecessor nodes in the graph (Vertex constrained by 2 lines, Vertex constrained by 1 face and 1 line, etc.). If a single view of a vertex is available, then labels are less restricted because the view ray can be used to help initialize the vertex. Given two views of a vertex, no restrictions are needed because even if the vertex is unconstrained, it can be initialized using view ray intersection (stereo). Similar restrictions are placed on line and face node labels.

These extra restrictions, which are referred to as "initiability" conditions, as used in step 40, basically insure that each element is completely specified. Satisfaction of these conditions guarantees that the system of equations being solved will not be singular due to some part of a model's structure being underspecified. Dependency graphs satisfying these conditions are also expected to produce better conditioned problems than graphs which do not, although this has not been proven.

Referring now to FIG. 5, a dependency graph 40' is illustrated for the building facade 24 of FIG. 2, as constrainted after reduction of the set of constraints specified in the "object house" constraint file, supra. Thus, a first graph row has input elements: the initial conditions specify the three Cartesian coordinate axes in graph elements 60-1 through 60-3; and the object frame origin in graph element 62. After these minimum spatial definitions are set, the other constraints are used in the second `line` of activity, to set the first vertex 27-1 at the origin (in graph element 64). Having now set this a priori condition for specifying lines, the second graph line of activity defines, in respective graph elements 66-1, 66-2 and 66-3, the directions of a first line (from point 1 at vertex 27-1 to point 2 at vertex 27-1) along the X axis, a second line (from point 1 to point 6 at vertex 27-6) along the Y axis, and a third line (from point 1 to point 4 at vertex 27-4) along the Z axis. The end points of the lines 12, 14 and 16 can now considered; the previous data can be used to define points 2 or 4, respectively at vertices 27-2 or 27-4, in graph element 68 or 70 of the fourth graph row. Because the specification of point 1 and lines 12 and 14 determine face 1234, just as face 6541 is determined by point 1 and lines 14 and 16, the fourth graph row also includes graph elements 72 and 74. Line 23 can now be found and placed in graph element 76 of row 5; in this particular dependency graph, row 5 also contains the graph elements 78 and 80 respectively for face and face 3478. This set of given and calculated constraints leads to the specification of line 43 (in element 82), line 4,7 (in element 84) and vertex 8 (in element 86) of the sixth graph row. Thereafter, the seventh row includes solutions for vertices 3 and 7 (in respective elements 88 and 90) and for line 54 (in element 92), while subsequent eighth, ninth and tenth (final) graphs rows respectively include respective elements 94, 96 and 98, which respectively specify vertex 5, line 65 and vertex 6. The entire model is now established for the object (building) having a selected point (vertex 1) included on its periphery.

The dependency graph is `written` by the graph search algorithm of my novel method. The inputs to the graph search algorithm are:

1. Scene topology (structure) graph--a graph structure provided by step 32 and including the object, face, line and vertex data structures, with links defining the connectivity between them.

2. List of constraints--each constraint (step 36 input) may define new point and direction vector geometric entities linked into the main scene structure graph.

3. Image feature location list--the measured positions in images of some but not necessarily all vertices, provided (at least in part) by the data input in step 34.

The outputs of the graph search algorithm are:

1. Dependency graph this step 40 data structure includes an assignment of one of the constrained geometric element types (from Table 1) to each of the geometric entities. Each element data structure is set up with address pointers to the elements upon which it depends.

2. Order list--a topological sort of the elements into a single list which preserves the partial ordering implied by the constraint dependencies. For an acyclic graph, a topological sort is always possible.

The graph search algorithm directly generates the order-listed topological sort of the elements by attempting to place the elements in order, one at a time. The first (presort) phase of the algorithm generates the topological sort in top-down manner, with least-dependent elements first. During the second phase, the final topological sort is generated in bottom-up, with most-dependent first. At any point during the search, the system seeks to add one more element to the order, and all elements will have been separated into two categories: those "determined" elements that appear after the next element to be placed in the graph; and those "undetermined" elements that appear before the next element to be placed.

During both phases of the algorithm, the system relies on "checkStatus" functions, associated with each of the types of geometric entities (face, line, vertex, frame, point, and direction vector). These functions are specialists at understanding the legal constraint relationships associated with each geometric entity. When a checkStatus function is invoked in the middle of the search, it is passed an address pointer to a particular instance of the associated geometric entity. The actions taken by these functions are:

1. Detect which connecting elements are constraining (according to the current state of the order list).

2. Ignore constraining elements whose constraints are subsumed by other elements. Example: a face's checkStatus function will ignore a constraining vertex if a line coincident with it is also constraining.

3. Decide if element is overconstrained.

4. Decide if element is underspecified (uninitiable).

5. Identify the element type if neither over- or underconstrained, by returning one of the type codes in the above Table.

The purpose of the initial top-down phase of the search is to establish global ordering information by determining the minimum depth in the dependency graph of each element. To accomplish this, a topological sort is generated, but elements are allowed to be overconstrained. The network is generated from top to bottom.

If any elements remain undetermined then stop and report underconstrained elements.

The purpose of the final bottom-up phase is to generate the final ordering of the elements, making sure that all constraints are compatible, and that no element is either over- or underconstrained. In this phase, the dependency graph, or network, is generated from the bottom up. The trial selection of elements is based on the depth number computed in the top-down phase.

After constraint elimination, the problem can be restated, in terms of the original nonlinear programming problem, as the placement of the parameters of the scene model, previously described by the k-dimensional vector x, into groups, where each group corresponds to a geometric entity of the scene (face, line, vertex, object orientation, etc.). A parameter group j is a function of the parameters of the elements on which it depends and on a (possibly null) set of free variables, y.sub.j.

Concatenating all the free variables y.sub.j, where j=1, . . . , n, creates a single free variable vector y, of dimension N equal to the number of degrees of freedom of the scene. All parameters of the geometric entities can be computed from the free variable vector x=x(y), and x is guaranteed to satisfy the constraint equations (1) and (2) for all y. The vertex positions are therefore functions of y, as is the function to be minimized. The optimization (1) and (2) is now equivalent to: ##EQU3## This optimization is unconstrained, and can be restated in somewhat more general terms by: concatenating all unknowns into one vector z; concatenating all modeled observations into one vector function g(z); concatenating all observation measurements into one vector g.sup.o ; and combining the standard deviations of observation errors into a diagonal weight matrix, W=diag(1/.sigma..sub.i). Then, the optimization becomes: ##EQU4##

The development of a robust procedure for use in step 40 to automatically generate an approximate solution to equation (4) is aided by the organization of the dependency graph. Considerable effort has gone into the development of a procedure to bootstrap from whatever a priori information is available (constraints plus measured vertex positions) to obtain initial estimates of unknown camera and scene parameters. No single procedure is capable of performing this initialization because it is highly dependent on what a priori or constraint information is available. Consequently, it is important to have several approaches, which are automatically drawn upon if they are applicable. The main difficulty is in estimating initial camera orientation. Algorithms are provided which could use any available control points, or direction constraints, or simply corresponding vertices between cameras to generate initial estimates that are usually accurate enough to allow convergence of the subsequent iteration.

The overall procedure for initializing the independent variables of the graph will be briefly summarized; once the camera models are approximated, the polygonal mesh is simple to initialize because each element is guaranteed to be initializable by the constraint elimination algorithm.

1. Unless already known, assume cameras located at infinity (c=0)

2. Initialize orientation, translation, and scale of cameras (if possible) using any "knowns" visible to the cameras:

lines in known directions

control points

3. Use stereo to approximate vertex locations (ignore constraints).

4. If some cameras are still unknown, go back to step 2.

5. Call each elements initialization procedure in an order such that all elements it depends on will be called first (always possible for an acyclic graph, which can be topologically sorted).

I have implemented step 42 by use of a Gauss-Newton iterative method for solving equation (4); this is the same as finding the least-squared-error solution to

These nonlinear equations are solved using the Gauss-Newton method, which iteratively improves an approximate solution z.sub.n. Each iteration of Gauss-Newton's method consists of:

Step 1: Compute Jacobian matrix, g'(z.sub.n)=.differential.g/.differential.z, which provides a linearized system of equations, substantially valid in the neighborhood of the z.sub.n.

Step 2: Solve the linearized system for the least-squared-error differential change, using the pseudoinverse (note that a simple adjustment to the g'.sup.T W.sup.2 g matrix in this step converts the Gauss-Newton algorithm to the Levenberg-Marquardt (LM) algorithm. The LM modification has a more robust convergence when used with a poor initial estimate, but in many cases takes significantly more steps to converge):

The pseudoinverse A.dagger. of matrix A is computed using:

Step 3: Apply the differential change to z.sub.n :

Actually, the sum found in step 3 only makes sense for certain z vector components, such as scalar camera parameters and point coordinates. Camera positions and much of the geometry, however, are represented by homogeneous transforms and to apply a differential change, [dp.sub.x dp.sub.y dp.sub.z d.phi..sub.x d.phi..sub.y d.phi..sub.z ].sup.T, to a transform, the change itself is converted into a homogeneous transform:

which then multiplies the original transform:

This technique for handling differential changes to a transform in the context of the Gauss-Newton procedure was previously described by Lowe.

Gauss-Newton's method converges to a local minimum mean squares error solution of Wg(z)=Wg.sup.o. The procedure described in previous section B provides an accurate enough starting point to allow convergence in most cases. For the times when the initial guess is not close enough (and convergence is not obtained), one may want to implement an additional solve mode, which is set to essentially use only the minimum model constraints needed to establish a base reference coordinate system, and simply computes camera models allowing the 3D points to float freely in space. A subsequent overall solve step then converges in all cases where a solution is well defined.

In my procedure, the computation of partial derivatives for the Jacobian begins with the camera model equations, which are derived hereinbelow in Appendix B. The camera model equations relate the image position vector i, to the vector position of a world point x.sup.w, depending on the value of the camera vector parameters f, R and p. ##EQU5## The point x.sup.c is the position of the point relative to the camera coordinate frame. A free parameter .lambda. of the model affects the sensor position through its effect on the variables x.sup.c and f. Hence, the evaluation of the partial derivatives of sensor coordinates with respect to the free parameter .lambda., will proceed along the lines of ##EQU6##

In this equation, x.sup.w refers to the world coordinates of the point, and C refers to the position of the camera coordinate frame. This equation (7) shows that the free parameter affects sensor position by its effect on the 3D camera coordinates of the point and by its effect on the internal camera parameters (focal length and scale in the initial implementation). Normally, a free parameter will affect either point coordinates or camera parameters and the other term drops out.

The two terms grouped in parenthesis in equation (7) indicate that the free parameter affects the camera coordinates of the point by either (1) its effect on the points position in the world, x.sup.w, or (2) its effect on the position of the camera with respect to the world. The evaluation of this last term requires concepts covered in Appendix C (relating to partial derivatives with respect to coordinate frames).

Several of the terms in equation (7) can be derived simply from equation (6): ##EQU7## Appendix C contains a derivation of the partial derivative of a point with respect to its base frame (which will be equation (12)): ##EQU8## The terms .differential.C/.differential..lambda. and .differential.f/.differential..lambda. evaluate to null vectors unless the free parameter is one of the camera parameters, in which case one element of one of the vectors is unity.

The one term in equation (7) that is not easily written is .differential.x.sup.w /.differential..lambda., giving the effect of the free parameter on the position of the vertex. Due to the chaining of elements in the dependency graph, a single free parameter can affect the positions of many vertices. In order to develop an automatic procedure for computing the overall Jacobian matrix, a small chain rule procedure is implemented for each of the 41 geometric elements in the above Table, which computes the partial derivative of the element with respect to an arbitrary free parameter in terms of the partial derivatives of all the terms it depends on. For example, for the face constrained by 1 vertex, the face is represented by a coordinate frame, F, and the vertex is represented by its world coordinate vector, v. Due to the constraint, ##EQU9## Given .differential.v/.differential..lambda., we can compute .differential.F/.differential..lambda. by following the procedure outlined in Appendix C: ##EQU10## Hence, using equation (8) in Appendix C, .delta.(.lambda.)=.theta. and d(.lambda.)=R.sup.T (v(.lambda.)-v(.lambda..sub.0)). Applying equation (13) in Appendix C: ##EQU11## This result is simple and easily written. Some of the other 40 cases are not so simple; however, they are all straight-forward to derive.

With these chain rule procedures in place for all geometric elements, the computation of the Jacobian matrix proceeds as follows:

In this procedure, in the loop where the partial derivatives of all affected elements are computed, the partial derivative functions are called in dependency order.

Referring now to FIGS. 6a-6d, several examples of objects extracted from reconstructed scenes are illustrated. The overall system used is PolyFit, a 3D model extraction system developed around the special-purpose nonlinear least-squares solver of the present invention, coupled with a user interface which allows a human user to employ the solver in an effective manner. Well over one hundred other models have so far been reconstructed, and the system is coming into routine use for building databases for the GE COMPU-SCENE.RTM. family of image generators. The imagery of the pillow-block 101 of FIG. 6a and the desk 102 scene of FIG. 6b was collected in my laboratory with a CCD video camera; two different views of the pillow block were taken, and eight different views of the desk scene were used. No additional information was used other than the images, and the initial constraints were those that a human operator easily surmised about the scene by examining the images. The building 103 of FIG. 6c was extracted from four images: one oblique view taken from a tower, a ground level handheld camera view, and two different overhead views. The building complex 104 of FIG. 6d was extracted from a total of 19 separate images. All of the views shown in FIGS. 6a-6d are hidden line wireframe renderings at viewpoints not corresponding to any of the images used to reconstruct the scene. In each of these examples, the scale of the models was specified by defining a length constraint within the scene. The base reference frame was arbitrarily attached to one of the objects in the scene. If a map which provides 3D positions of recognizable points were available, then the models could be reconstructed and scaled within a coordinate frame related to the map. These examples illustrate the wide variety of models that can be reconstructed with the system and method of the present invention.

My novel system and methods for extracting object data from scene geometry reconstructed from images, using a maximum-likelihood estimation to account for random measurement errors, has been described herein in great detail. Specifically, the maximum-likelihood scene was shown to be the solution to a constrained nonlinear least square minimization, for which a solution procedure was developed. Many modifications and variation will now occur to those skilled in the art. It is my intent to be limited only by the scope of the appending claims and not by the details and instrumentalities presented herein by way of explanation of the particular embodiment presently preferred.

In the implementation of this system, extensive use is made of the homogeneous transform (rotation matrix and translation vector) which represents the position of one certainty frame with respect to another. In addition to being the representation of various geometric elements (planes, lines, directions), homogeneous transforms are used to describe the scene geometry as a two-level hierarchy, with the position of each object in the world defined by a homogeneous transform, and the internal geometry of each object defined within its own certainty frame. One reason that each object has its own coordinate frame is that constraints on the internal geometry of objects are most conveniently defined relative to a preferred-object-oriented coordinate frame. Also, having an object reference frame for each object effectively separates the internal geometry of the object from the orientation and position of the object within the scene. The system can then, for instance, take fixed object models and position them within the scene.

The following are useful ways of referring to the components of the 4.times.4 matrix homogeneous transform T: ##EQU12##

Camera modeling schemes normally place the camera coordinate frame at either the optical focal point or the focal plane. This choice leads to two potential problems for iterative searches (which problems have been alleviated with the novel camera modeling system and methods of the present invention). The first problem arises if the camera coordinate frame is placed in proximity to where the camera is physically located. The problem is (assuming for the sake of discussion that the Z-axis is aligned with the optical axis) that the first order effect of rotations about the X-axis and translations along the Y-axis are the same: the image translates in the image Y direction. This first order equivalence is also true of rotations about the Y-axis and translations along the X-axis. The second problem is that the first-order effect of a change in focal length and a change in camera range (translation along the Z-axis) are the same whenever the camera range is many times greater than the depth of the scene; the perspective transformation reduces to orthographic projection, and the effect of both range and focal length are a change in scale. The effect of these ambiguities in rotations/translations and range/focal length is to increase the condition number of the Jacobian matrix, and gradient-based methods will hem-stitch. Each step of the iteration must be chosen to balance the large first-order effects, hence the system becomes less sensitive to the second-order effects (perspective and rotation of the points within the scene).

The first problem is avoided by choosing to represent the camera position with a coordinate frame called the local camera frame, which lies on the optical axis of the camera optics, but which is located in the proximity of the scene. This local frame will be oriented the same as the remote frame which lies at the focal center of the optics. The Z-axis of both frames is directed away from the scene (hence all points in the scene will have negative z values in the remote frame). The X-axis is directed in the direction of increasing raster lines and the Y-axis is directed in the direction of a left-to-right scan of the image. Note that rotations of this axis will have the desired effect of rotating the camera about the scene while the scene remains in view of the camera.

The second problem is avoided by formulating the perspective transformation such that the two parameters "range" and "focal length" are replaced with two new parameters "closeness" (inverse range) and "scale" (ratio of focal length to range).

As shown in FIG. 7, in a camera model which avoids both the first (remote rotations) problem and the second (special orthographic projection case) problem, a point (x.sup.r, y.sup.r, z.sup.r) in the remote camera frame (located at the camera focal center) is projected to the point (u,v) in the focal plane (located at a distance f from the focal center), according to the perspective transformation:

To avoid the first problem in this form of the camera model, I define the position of the point with respect to a local camera frame (x.sup.c,y.sup.c,z.sup.c) where (x.sup.c,y.sup.c)=(x.sup.r,y.sup.r),z.sup.c =z.sup.r +1/c. The second problem is avoided by the use of a parameter `c` to represent the inverse range from the local frame to the remote frame along the optical axis. With these substitutions, the perspective transformation becomes:

Lastly, define a scale parameter, s=fc, which scales distances from the local frame to the focal plane (in the absence of perspective distortion). The final result is:

Note that this formulation, based on closeness and scale parameters, reduces cleanly to an orthographic projection as the camera range goes to infinity, because the closeness parameter `c` goes to a well-defined point (zero). However, note also that now there is an extra parameter, due to the fact that the position of the local coordinate frame along the optical axis is arbitrary. To handle the ambiguity associated with the extra parameter, the position of the local frame along the optical axis will be held fixed at each step of the iteration, and the position of the remote frame moved along the axis by varying the closeness parameter.

The camera model also normally includes four parameters (a.sub.u, b.sub.u, a.sub.v, b.sub.v) that relate pixel position g,h in the digitized image to the u,v position in the focal plane. The relevant equations are:

The `a` factors are sensor scale factors, in units of pixels/unit distance, and the `b` factors are the principal point offset in pixels.

Summarizing, the camera parameters c include a homogeneous transform (rotation matrix R and translation p) representing the position of the word with respect to the camera local frame, closeness c, scale s, and internal sensor parameters a.sub.u, b.sub.u, a.sub.v, b.sub.v. The camera projection function p introduced earlier is defined by: ##EQU13##

1. Partial Derivative with respect to (wrt) Coordinate Frames:

To obtain partial derivatives with respect to rotations and translations of coordinate frames, I work with differential transforms (see, e.g. Robot Manipulators by Richard Paul). A differential transform consists of a differential translation, represented by a vector d=[d.sub.x d.sub.y d.sub.z ].sup.T, and a differential rotation, represented by a vector .delta.=[.delta..sub.x .delta..sub.y .delta..sub.z ].sup.T. A differential transform can be represented by stacking the translation and rotation vectors into a single 6-vector:

In practices, a homogeneous transform Diff is obtained from these vectors by:

Here, Trans(d.sub.x,d.sub.y,d.sub.z) is a pure translation transform, and Rot(.delta..sub.x,.delta..sub.y,.delta..sub.z) is a rotation about the X-axis by .delta..sub.x, followed by a rotation about the Y-axis by .delta..sub.y, followed by a rotation about the Z-axis by .delta..sub.z.

For small values of d and .delta., the result of applying the Diff transform to a vector x, can be approximated by using the vector cross product to perform the rotation: ##EQU14## Hence Diff can also be expressed as ##EQU15## where the notation [.delta.].sub.x defines a mapping from a 3D vector to a 3-by-3 matrix: ##EQU16## Using this notation, the cross production of two vectors can be expressed as multiplication of a 3-by-3 matrix and a vector:

Which expression for the differential transform is most convenient to use depends on the task at hand. If I am applying the differential transform computationally, I use the homogeneous transform expression (8), because it will preserve the orthogonality of rotation matrices. If I am analytically evaluating differentials, the cross product form (9) will usually be the most useful.

The partial derivative of some function of a coordinate frame T, with respect to changes in that frame, can be found using: ##EQU17## 2. Partial Derivative of Transformed Vector wrt the Transformation:

As an example of use of equation (11), suppose I have y=f(T)=Tx=Rx+d. It should be noted that, in order to avoid the use of homogeneous vectors, the homogeneous transforms such as T and Diff, although written as 4-by-4 matrices, are treated as operators that transform 3D-vectors into 3D-vectors by a rotation and translation. Then, making use of equation (9): ##EQU18## I now apply equation (11): ##EQU19## In equation (12), x, y and z are the coordinates of x (with respect to T). 3. Partial Derivative of a Point wrt its Base Frame:

The result in equation (12) suggested to me a definition for the partial derivative x of a point, with respect to a differential change in its base transformation T, given by: ##EQU20## With this definition of .differential.x/.differential.T, equation (12) can be derived using the chain rule: ##EQU21## while recognizing that .differential.y/.differential.x=R, and expanding .differential.x/.differential.T using equation (13).

4. Partial Derivative of Coordinate Frame wrt Free Parameter:

I assume a homogeneous transform T(.lambda.), which varies as a function of a free parameter .lambda.. To compute the partial derivative of the transform with respect to the free parameter, I expressed the changes in T in terms of a differential transform. Thus, in the neighborhood of .lambda.=.lambda..sub.0, the variation of T can be represented by a differential transform Diff(.lambda.), as follows:

As defined in equation (8), the differential transform Diff depends on the 6-vector D. The partial derivative of the coordinate frame (ie, the transform) with respect to .lambda. is then defined to be: ##EQU22## In practice, the partial derivatives will often be determined by simply considering the translational and rotational effects of small changes in the free parameter. In other contexts, the partial derivatives will be determined using the chain rule.

5. Partial Derivative of Coordinate Frame wrt Coordinate Frame:

To calculate the partial derivative of one coordinate frame with respect to a second coordinate frame which is treated as a free parameter, I express the differential changes induced in the first frame by differential changes in the second frame. I let the functional relationship be T(A). First, I establish the relationship between differential changes in A on T:

And now the partial derivative is defined by: ##EQU23## 6. Transforming Differential Transforms:

At least one formula is known for transforming a differential transform from one coordinate frame to another. The version of this formula which I have used is derived using equations (15) and (16). To do this, I compute the partial derivative of T=A C with respect to A. From equation (15), I obtain:

From the above, I obtain the following expression for Diff.sub.T :

Expanding, letting C be written in terms of a rotation matrix R and translation vector p, gives: ##EQU24## Using equation (9), I can equate this expression with Diff.sub.T in terms of its differential parameters: ##EQU25## From the upper left partition of this matrix equality, we obtain: ##EQU26## The last step required the formula:

which is true for any orthogonal matrix R. This formula can be verified by carrying out the multiplication of the three terms and simplifying. The truth of the formula can be simply recognized by considering each terms as a transformation. The transformation on the left hand side consists of a rotation followed by taking the cross product with a vector .delta., and then rotating the result back to the original coordinate frame. These operations are equivalent to taking the cross product with the vector .delta. rotated into the original coordinate frame, which is the transformation on the right hand side.

From the upper right partition of matrix equality (17): ##EQU27## The notation p.times.R.tbd.[p].sub.x R (a vector crossed with a matrix) represents the matrix with each column equal to p crossed with the corresponding column of R. Combining equations (18) and (19), we obtain an expression for the transformation of differential changes in A into differential changes in T=AC: ##EQU28## Finally, applying equation (16), we obtain the partial derivative of T with respect to A: ##EQU29##