Stereo reconstruction employing a layered approach and layer refinement techniques

A system and method for extracting structure from stereo that represents the scene as a collection of planar layers. Each layer optimally has an explicit 3D plane equation, a colored image with per-pixel opacity, and a per-pixel depth value relative to the plane. Initial estimates of the layers are made and then refined using a re-synthesis step which takes into account both occlusions and mixed pixels. Reasoning about these effects allows the recovery of depth and color information with high accuracy, even in partially occluded regions. Moreover, the combination of a global model (the plane) with a local correction to it (the per-pixel relative depth value) imposes enough local consistency to allow the recovery of shape in both textured and untextured regions.

BACKGROUND
 1. Technical Field
 The invention is related to a system and process for extracting 3D
 structure from plural, stereo, 2D images of a scene by representing the
 scene as a group of image layers characterized by estimated parameters
 including the layer's orientation and position, per-pixel color, per-pixel
 opacity, and optionally a residual depth map, and more particularly, to
 such a system and process for refining the estimates for these layer
 parameters.
 2. Background Art
 Extracting structure from stereo has long been an active area of research
 in the imaging field. However, the recovery of pixel-accurate depth and
 color information from multiple images still remains largely unsolved.
 Additionally, existing stereo algorithms work well when matching feature
 points or the interiors of textured objects. However, most techniques are
 not sufficiently robust and perform poorly around occlusion boundaries and
 in untextured regions.
 For example, a common theme in recent attempts to solve these problems has
 been the explicit modeling of the 3D volume of the scene. The volume of
 the scene is first discretized, usually in terms of equal increments of
 disparity. The goal is then to find the so-called voxels which lie on the
 surfaces of the objects in the scene using a stereo algorithm. The
 potential benefits of these approaches can include, the equal and
 efficient treatment of a large number of images, the explicit modeling of
 occluded regions, and the modeling of mixed pixels at occlusion boundaries
 to obtain sub-pixel accuracy. However, discretizing space volumetrically
 introduces a huge number of degrees of freedom. Moreover, modeling
 surfaces by a discrete collection of voxels can lead to sampling and
 aliasing artifacts.
 Another active area of research directed toward solving the aforementioned
 problems is the detection of multiple parametric motion transformations
 within image sequence data. The overall goal is the decomposition of the
 images into sub-images (or "layers") such that the pixels within each
 sub-image move consistently with a single parametric transformation.
 Different sub-images are characterized by different sets of parameter
 values for the transformation. A transformation of particular importance
 is the 8-parameter homography (collineation), because it describes the
 motion of points on a rigid planar patch as either it or the camera moves.
 The 8 parameters of the homography are functions of the plane equations
 and camera matrices describing the motion.
 While existing layer extraction techniques have been successful in
 detecting multiple independent motions, the same cannot be said for scene
 modeling. For instance, the fact that the plane equations are constant in
 a static scene (or a scene imaged by several cameras simultaneously) has
 not been exploited. This is a consequence of the fact that, for the most
 part, existing approaches have focused on the two frame problem. Even when
 multiple frames have been considered, it has primarily been solely for the
 purposes of using past segmentation data to initialize future frames.
 Another important omission is the proper treatment of transparency. With a
 few exceptions, the decomposition of an image into layers that are
 partially transparent (translucent) has not been attempted.
 SUMMARY
 The present invention relates to stereo reconstructions that optimally
 recover pixel-accurate depth and color information from multiple images,
 including around occlusion boundaries and in untextured regions. This is
 generally accomplished using an approach to the stereo reconstruction that
 represents the 3D scene as a collection of approximately planar layers,
 where each layer has an explicit 3D plane equation and a layer sprite
 image, and may also be characterized by a residual depth map. The layer
 sprite refers to a colored image with a defined per-pixel opacity
 (transparency). The residual depth map refers to a per-pixel depth value
 relative to the plane. The approach of segregating the scene into planar
 components allows a modeling of a wider range of scenes. To recover the
 structure of the scene, standard techniques from parametric motion
 estimation, image alignment, and mosaicing can be employed.
 More specifically, the approach to the stereo reconstruction based on
 representing the 3D scene as a collection of approximately planar layers
 involves first estimating the desired parameters (e.g. plane equation,
 sprite image and depth map) and then refining these estimates. The
 estimating phase can be accomplished via any appropriate method, such as
 the methods disclosed in a co-pending application entitled STEREO
 RECONSTRUCTION EMPLOYING A LAYERED APPROACH by the inventors of this
 application and assigned to the common assignee. This application was
 filed on Mar. 20, 1998 and assigned Ser. No. 09/045,519. The full approach
 disclosed in the application, which is believed to provide the best
 estimate of the layer parameters, and so the structure of the 3D scene,
 includes:
 (a) inputting plural 2D images as well as camera projection matrices
 defining the location and orientation of the camera(s) responsible for
 creating each image, respectively;
 (b) assigning each pixel making up each 2D image to one of the plural
 layers;
 (c) estimating a plane equation for each layer that defines the orientation
 and position of that layer in 3D space;
 (d) estimating a sprite image for each layer characterized by a per-pixel
 color and a per-pixel opacity;
 (e) estimating a residual depth map for each layer wherein each residual
 depth map defines the distance each pixel of the associated layer is
 offset from the estimated plane of that layer;
 (f) re-estimating each layer's sprite image based on the residual depth map
 associated with the layer;
 (g) re-assigning pixels assigned to a particular layer to another layer by
 using the estimates for the plane equation, sprite image, and residual
 depth map for each layer as a guide;
 (h) iteratively repeating steps (c) through (g) for each layer until the
 change in the value of at least one layer parameter relative to its value
 in an immediately preceding iteration falls below a prescribed threshold
 assigned to the parameter; and
 (i) outputting data representative of the plane equation, sprite image and
 residual depth map estimates for each layer.
 Only the input, pixel assignment, plane equation and sprite image
 estimation, and output modules (less the residual depth map) are necessary
 to produce a useable layered representation of the scene. However, the
 accuracy of the layered representation can be progressively improved with
 the respective addition of each of the remaining modules, i.e. the depth
 map estimation, sprite image re-estimation, and pixel re-assignment and
 iteration modules.
 The initial estimates of the layer parameters are refined in accordance
 with systems and methods embodying the present invention. The refining
 process is accomplished using a re-synthesis, which takes into account
 both occlusions and mixed pixels. Approximate knowledge of the 3D
 structure (camera matrices and plane equations) allows reasoning about the
 image formation process. Specifically, a forward (generative) model of
 image synthesis is used, as well as a process of measuring of how well the
 layers re-synthesize the input images. Optimizing this measure allows
 refinement of the layer sprite estimates, and, in particular, the estimate
 of their true colors and opacities. This approach results in the correct
 recovery of mixed pixels, a step which is necessary to obtain sub-pixel
 accuracy and to ensure robustness at occlusion boundaries. Once the layer
 sprite estimates are refined, the plane equations and residual depth maps
 (if employed) can be refined as well using the original estimation
 process, such as the one disclosed in the aforementioned co-pending
 application.
 The layered approach to stereo reconstruction shares many of the advantages
 of the previously described volumetric approaches because the 3D
 information contained in the layers is used to reason about occlusion and
 mixed pixels. However, the layered approach according to the present
 invention offers a number of additional advantages, including:
 A combination of the global model (the plane) and the local correction to
 it (the per-pixel depth map) that results in very robust performance and
 extremely accurate depth maps;
 A layered approach that enables the recovery of scene structure in
 untextured regions because there is an implicit assumption that untextured
 regions are planar--this is not an unreasonable assumption, especially in
 man-made environments;
 A form of the output (i.e., a collection of approximately planar regions
 with per-pixel depth offsets) that is more suitable than a discrete
 collection of voxels for many applications, including, view interpolation
 and interactive scene modeling for rendering and video parsing.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS
 In the following description of the preferred embodiments of the present
 invention, reference is made to the accompanying drawings which form a
 part hereof, and in which is shown by way of illustration specific
 embodiments in which the invention may be practiced. It is understood that
 other embodiments may be utilized and structural changes may be made
 without departing from the scope of the present invention.
 FIG. 1 and the following discussion are intended to provide a brief,
 general description of a suitable computing environment in which the
 invention may be implemented. Although not required, the invention will be
 described in the general context of computer-executable instructions, such
 as program modules, being executed by a personal computer. Generally,
 program modules include routines, programs, objects, components, data
 structures, etc. that perform particular tasks or implement particular
 abstract data types. Moreover, those skilled in the art will appreciate
 that the invention may be practiced with other computer system
 configurations, including hand-held devices, multiprocessor systems,
 microprocessor-based or programmable consumer electronics, network PCs,
 minicomputers, mainframe computers, and the like. The invention may also
 be practiced in distributed computing environments where tasks are
 performed by remote processing devices that are linked through a
 communications network. In a distributed computing environment, program
 modules may be located in both local and remote memory storage devices.
 With reference to FIG. 1, an exemplary system for implementing the
 invention includes a general purpose computing device in the form of a
 conventional personal computer 20, including a processing unit 21, a
 system memory 22, and a system bus 23 that couples various system
 components including the system memory to the processing unit 21. The
 system bus 23 may be any of several types of bus structures including a
 memory bus or memory controller, a peripheral bus, and a local bus using
 any of a variety of bus architectures. The system memory includes read
 only memory (ROM) 24 and random access memory (RAM) 25. A basic
 input/output system 26 (BIOS), containing the basic routine that helps to
 transfer information between elements within the personal computer 20,
 such as during start-up, is stored in ROM 24. The personal computer 20
 further includes a hard disk drive 27 for reading from and writing to a
 hard disk, not shown, a magnetic disk drive 28 for reading from or writing
 to a removable magnetic disk 29, and an optical disk drive 30 for reading
 from or writing to a removable optical disk 31 such as a CD ROM or other
 optical media. The hard disk drive 27, magnetic disk drive 28, and optical
 disk drive 30 are connected to the system bus 23 by a hard disk drive
 interface 32, a magnetic disk drive interface 33, and an optical drive
 interface 34, respectively. The drives and their associated
 computer-readable media provide nonvolatile storage of computer readable
 instructions, data structures, program modules and other data for the
 personal computer 20. Although the exemplary environment described herein
 employs a hard disk, a removable magnetic disk 29 and a removable optical
 disk 31, it should be appreciated by those skilled in the art that other
 types of computer readable media which can store data that is accessible
 by a computer, such as magnetic cassettes, flash memory cards, digital
 video disks, Bernoulli cartridges, random access memories (RAMs), read
 only memories (ROMs), and the like, may also be used in the exemplary
 operating environment.
 A number of program modules may be stored on the hard disk, magnetic disk
 29, optical disk 31, ROM 24 or RAM 25, including an operating system 35,
 one or more application programs 36, other program modules 37, and program
 data 38. A user may enter commands and information into the personal
 computer 20 through input devices such as a keyboard 40 and pointing
 device 42. Of particular significance to the present invention, a camera
 55 (such as a digital/electronic still or video camera, or
 film/photographic scanner) capable of capturing a sequence of images 56
 can also be included as an input device to the personal computer 20. The
 images 56 are input into the computer 20 via an appropriate camera
 interface 57. This interface 57 is connected to the system bus 23, thereby
 allowing the images to be routed to and stored in the RAM 25, or one of
 the other data storage devices associated with the computer 20. However,
 it is noted that image data can be input into the computer 20 from any of
 the aforementioned computer-readable media as well, without requiring the
 use of the camera 55. Other input devices (not shown) may include a
 microphone, joystick, game pad, satellite dish, scanner, or the like.
 These and other input devices are often connected to the processing unit
 21 through a serial port interface 46 that is coupled to the system bus,
 but may be connected by other interfaces, such as a parallel port, game
 port or a universal serial bus (USB). A monitor 47 or other type of
 display device is also connected to the system bus 23 via an interface,
 such as a video adapter 48. In addition to the monitor, personal computers
 typically include other peripheral output devices (not shown), such as
 speakers and printers.
 The personal computer 20 may operate in a networked environment using
 logical connections to one or more remote computers, such as a remote
 computer 49. The remote computer 49 may be another personal computer, a
 server, a router, a network PC, a peer device or other common network
 node, and typically includes many or all of the elements described above
 relative to the personal computer 20, although only a memory storage
 device 50 has been illustrated in FIG. 1. The logical connections depicted
 in FIG. 1 include a local area network (LAN) 51 and a wide area network
 (WAN) 52. Such networking environments are commonplace in offices,
 enterprise-wide computer networks, intranets and the Internet.
 When used in a LAN networking environment, the personal computer 20 is
 connected to the local network 51 through a network interface or adapter
 53. When used in a WAN networking environment, the personal computer 20
 typically includes a modem 54 or other means for establishing
 communications over the wide area network 52, such as the Internet. The
 modem 54, which may be internal or external, is connected to the system
 bus 23 via the serial port interface 46. In a networked environment,
 program modules depicted relative to the personal computer 20, or portions
 thereof, may be stored in the remote memory storage device. It will be
 appreciated that the network connections shown are exemplary and other
 means of establishing a communications link between the computers may be
 used.
 The exemplary operating environment having now been discussed, the
 remaining part of this description section will be devoted to a
 description of the program modules embodying the invention and the testing
 of these modules. Unless stated otherwise, the following description
 employs homogeneous coordinates for both 3D world coordinates x=(x, y, z,
 1).sup.T and for 2D image coordinates u=(u, v, 1).sup.T. In addition, it
 is assumed that when indexing 2D images, the homogeneous coordinates are
 converted into 2D pixel indices via conventional methods before accessing
 the array containing the image. Referring now to FIG. 2, suppose K images
 I.sub.1 (u.sub.1), I.sub.2 (u.sub.2), . . . , I.sub.K (u.sub.K) captured
 by K cameras with projection matrices P.sub.1, P.sub.2, . . . , P.sub.K.
 In what follows, the image coordinates u.sub.k will be dropped unless they
 are needed to explain a warping operation explicitly.
 The basic concept is that the world can be reconstructed as a collection of
 L approximately planar layers. A sprite image of a layer is denoted by
 L.sub.l (u.sub.l)=(.alpha..sub.l.sup..multidot. r.sub.l,
 .alpha..sub.l.sup..multidot. g.sub.l, .alpha..sub.l.sup..multidot.
 b.sub.l, .alpha..sub.l) where r.sub.l =r.sub.l (u.sub.l) is the red band,
 g.sub.l =g.sub.l (u.sub.l) is the green band, b.sub.l =b.sub.l (u.sub.l)
 is the blue band, and .alpha..sub.l =.alpha..sub.l (u.sub.l) is the
 opacity of pixel u.sub.l. Each layer is also associated with a homogeneous
 vector n.sub.l, which defines the plane equation of the layer via
 n.sub.l.sup.T X=0, and optionally a per-pixel residual depth offset
 Z.sub.l (u.sub.l). Thus, the scene is represented by L sprite images
 L.sub.l on planes n.sub.l.sup.T =0 with depth offsets Z.sub.l.
 Techniques for estimating the layer parameters n.sub.l, L.sub.l, and
 Z.sub.l, typically require an assumption of boolean valued opacities where
 a pixel in an image is characterized as being either completely opaque or
 completely transparent. The aforementioned co-pending application is an
 example where such an assumption is made for estimation purposes. However,
 the techniques embodying the present invention use these estimates of the
 structure as an input into a re-synthesis process that refines the layer
 sprites L.sub.l, to including real opacity values .alpha..sub.l.
 By way of an example, the methods employed in the previously-identified
 co-pending application will be described herein to show how the layer
 parameters L.sub.l, n.sub.l, and Z.sub.l are first estimated prior to
 employing the refinement techniques associated with the present invention.
 The estimation of the layer parameters can be subdivided into a number of
 steps for respectively estimating L.sub.l, n.sub.l, and Z.sub.l. To
 compute these quantities, we need to add auxiliary boolean mask images
 B.sub.kl. The boolean masks B.sub.kl denote the pixels in image I.sub.k
 which are images of points in layer L.sub.l. In other words, the mask
 images B.sub.kl define which pixels are in a particular layer. Since we
 are assuming boolean opacities, B.sub.kl =1 if and only if L.sub.l is the
 front-most layer which is opaque at that pixel in image I.sub.k. Hence, in
 addition to L.sub.l, n.sub.l, and Z.sub.l, we also need to estimate the
 boolean masks B.sub.kl. Once we have estimated these masks, we can compute
 masked input images M.sub.kl =B.sub.kl.TM.I.sub.k.
 FIG. 3 illustrates a high-level structure of a preferred layered stereo
 reconstruction approach employing refinement of the estimated layer
 parameters L.sub.l, n.sub.l, and Z.sub.l, using real opacity values
 .alpha..sub.l. The estimation steps associated with the aforementioned
 co-pending application have been included as the preferred method for
 estimating the layer parameters L.sub.l, n.sub.l, and Z.sub.l. However,
 other layer estimation techniques could also be employed if desired.
 Referring to FIG. 3, the input consists of the aforementioned collection
 of images I.sub.k taken with cameras of known geometry P.sub.k. The goal
 is to estimate the layer sprites L.sub.l, the plane vectors n.sub.l, and
 preferably the residual depths Z.sub.l. Given any three of L.sub.l,
 n.sub.l, Z.sub.l, and B.sub.kl, there are techniques for estimating the
 remaining one. The approach generally consists of first initializing these
 quantities. Then, each of these quantities is in turn iteratively
 estimated while fixing the other three. After good initial estimates of
 the layers are obtained, the refinement process embodying the present
 invention is begun in which real valued opacities .alpha..sub.l are used
 and the layer sprites are refined, including the opacities. In the
 following discussion, each of the steps in the layer parameter estimation
 process of the previously-identified co-pending application is discussed.
 Then, the subject process involving an image formation model and a
 re-synthesis procedure used to refine each of the layer sprites is
 discussed.
 Additionally, in the remainder of this specification, the description
 refers to various individual publications identified by an alphanumeric
 designator contained within a pair of brackets. For example, such a
 reference may be identified by reciting, "reference [WA93]" or simply
 "[WA93]". Multiple references will be identified by a pair of brackets
 containing more than one designator, for example, [WA93, SA96, BJ96]. A
 listing of the publications corresponding to each designator can be found
 at the end of this Detailed Description section.
 1. Computation of Layers Using Boolean Opacities
 1.1 Initialization of Layers
 Initialization of the layers is a process by which the pixels forming the
 images are initially assigned to a particular layer, thereby segmenting
 the images. The Boolean mask images B.sub.kl, and so the masked input
 images M.sub.kl are defined from this initial assignment of pixels. There
 are a number of methods of initializing the layers. One way is to
 initialize a large number of small random layers, which are then allowed
 to grow and merge until a small number of larger layers remain which
 accurately model the scene. Such an approach was employed in references
 [WA93, SA96, BJ96]. A second approach would be to use a sequential method
 in which a dominant motion is first extracted and then the residual
 regions are processed recursively. This method was employed in references
 [IAH95, SA96]. A third possibility would be to perform a color or
 intensity segmentation in each image as used in reference [ASB94]. The
 resulting segments could be matched and used as initial layer assignments.
 A fourth approach consists of first applying a conventional stereo
 algorithm to get an approximate depth map. A plane fitting algorithm could
 then be used to initialize the layers. A final method would be to get a
 user to initialize the layers manually. In many applications, such as
 model acquisition [DTM96] and video parsing [SA96], the goal is a
 semi-automatic process where limited user input is acceptable.
 1.2 Estimation of Plane Equations
 Referring to FIG. 4, in order to compute the plane equation vector n.sub.l,
 it is necessary to map points in the masked images M.sub.kl, derived as
 previously discussed (step 400), onto the plane n.sub.l.sup.T X=0. If X is
 a 3D world coordinate of a point and u.sub.k is the image of X in camera
 P.sub.k :
EQU u.sub.k =P.sub.k X (1)
 where equality is in the 2D projective space .rho..sup.2. Since P.sub.k is
 of rank 3, it can be written:
EQU X=P.sub.k *u.sub.k +sp.sub.k (2)
 where P.sub.k *=P.sup.T.sub.k (P.sub.k p.sup.T.sub.k).sup.-1 is the
 pseudoinverse of P.sub.k, s is an unknown scalar, and p.sub.k is a vector
 in the null space of P.sub.k, i.e. P.sub.k p.sub.k =0. If X lies on the
 plane n.sub.l.sup.T X=0 then:
 n.sub.l.sup.T P.sub.k *u.sub.k +sn.sub.l.sup.T p.sub.k =0 (3)
 Solving this equation for s, substituting into Equation (2), and
 rearranging using the fact that equality is up to a scale factor in
 .rho..sup.2, yields:
EQU X=((n.sub.l.sup.T p.sub.k)I-p.sub.k n.sub.l.sup.T)P.sub.k *u.sub.k (4)
 The importance of Equation (4) is that it facilitates the mapping of a
 point u.sub.k in image M.sub.kl onto the point on plane n.sub.l.sup.T X=0,
 of which it is an image. Afterwards this point is mapped onto its image in
 another camera P.sub.k' :
EQU u.sub.k' =P.sub.k' ((n.sub.l.sup.T p.sub.k)P.sub.k
 *u.sub.k.ident.H.sup.l.sub.kk' u.sub.k (5)
 where H.sup.l.sub.kk ' is a homography (collineation of .rho..sup.2)
 Equation (5) describes the image coordinate warp between the two images
 M.sub.kl and M.sub.k'l, which would hold if all the masked image pixels
 were images of world points on the plane n.sub.l.sup.T X 0. Using this
 relation, all of the masked images can be warped onto the coordinate frame
 of one distinguished image, M.sub.1l (step 402), without loss of
 generality, as follows:
EQU (H.sup.l.sub.1 k o M.sub.kl)(u.sub.i).ident.M.sub.kl (H.sup.l.sub.1 k
 u.sub.1) (6)
 Here, H.sup.l.sub.1 k o M.sub.kl is the masked image M.sub.kl, warped into
 the coordinate frame of M.sub.1l.
 It should be noted that it would also be possible to add an extra 2D
 perspective coordinate transformation in the just-described method (not
 shown in FIG. 4). Suppose H is an arbitrary homography. Each masked image
 could be warped onto H o H.sup.l.sub.1 k o M.sub.kl
 (u.sub.1).ident.M.sub.kl (HH.sup.l.sub.1k U.sub.1). The addition of the
 homography H can be used to remove the dependence of the above on one
 distinguished image. For example, H could be chosen to optimize some
 function of the homographies HH.sup.l.sub.1 k which measures how close
 they are to identity.
 The property used to compute n.sub.l is the following: assuming the pixel
 assignments to layers B.sub.kl are correct, the world is piecewise planar,
 and the surfaces are Lambertian, the warped images H.sup.l.sub.1 k o
 M.sub.kl should agree with each other and with M.sub.1l where they
 overlap. There are a number of functions which could possibly be used to
 measure the degree of consistency between the warped images (step 404),
 including least squares measures [BAHH92] and robust measures [DP95,
 SA96]. In both cases, the goal is the same, find the plane equation vector
 n.sub.l which maximizes the degree of consistency (step 406), which is
 usually a numerical minimum. Typically, this extremum is found using some
 form of gradient decent, such as the Gauss-Newton method, and the
 optimization is performed in a hierarchical (i.e. pyramid based) fashion
 to avoid local extrema. To apply this standard approach, simply derive the
 Jacobian of the image warp H.sup.l.sub.1 k with respect to the parameters
 of n.sub.l as described in reference [SS98]. This can be computed in a
 straightforward manner from Equation (5) because the camera matrices
 P.sub.k are assumed to be known. It is noted that the plane defined by
 n.sub.l can turn out to be of any orientation and need not be a frontal
 plane.
 1.3 Estimation of Layer Sprites Images
 Referring to FIG. 5, before the layer sprite images L.sub.l can be
 computed, 2D coordinate systems must be chosen for the planes. Such
 coordinate systems can be specified by a collection of arbitrary (rank 3)
 camera matrices Q.sub.l. An interesting choice for Q.sub.l is one for
 which the null space of Q.sub.l is perpendicular to the plane, and for
 which the pseudo-inverse maps the coordinate axes onto perpendicular
 vectors in the plane (i.e., a camera for which the plane is a frontal
 image plane). However, it is often undesirable to use a fronto-parallel
 camera since this may unnecessarily warp the input images.
 Once Q.sub.l has been specified, it can be shown, by the same argument used
 in Equations (4) and (5), that the image coordinates u.sub.k of the point
 in image M.sub.kl which is projected onto the point u.sub.l on the plane
 n.sub.l.sup.T X=0 is given by:
EQU u.sub.k =P.sub.k ((n.sub.l.sup.T q.sub.l)I-q.sub.l n.sub.l.sup.T)Q.sub.l
 *u.sub.l.ident.H.sup.l.sub.k u.sub.k (7)
 where Q.sub.l * is the pseudo-inverse of Q.sub.l and q.sub.l is a vector in
 the null space of Q.sub.l. The homography H.sup.l.sub.k warps the
 coordinate frame of the plane backward onto that of image M.sub.kl. The
 homography can also be used to warp the image M.sub.kl forward onto the
 plane (step 500), the result of which is denoted H.sub.l.sub.k o M.sub.kl.
 After the masked image has been warped onto the plane, the layer sprite
 image (with boolean opacities) can be estimated by "blending" the warped
 images:
 ##EQU1##
 where .sym. is the blending operator (step 502).
 There are a number of ways the blending could be performed. One simple
 method would be to take the mean of the color or intensity values. A
 refinement would be to use a "feathering" algorithm, where the averaging
 is weighted by the distance of each pixel from the nearest invisible pixel
 in M.sub.kl, as used in reference [SS97]. Alternately, robust techniques
 could be used to estimate L.sub.l from the warped images. The simplest
 example of such a technique is the median, but many others exist.
 There is a possibility that the images I.sub.k could have different
 exposure characteristics owing to differences in the cameras used to
 capture the images or differences in the camera settings between images if
 the same camera was used to take the disparate images. The difference in
 the exposure between images will skew the colors and opacities exhibited
 by corresponding image points. One possibility for addressing this
 exposure problem is to simply ignore it. The previously-described blending
 operation will tend to smooth out the disparities between the images. In
 some applications, this smoothing will produce sprite images reasonably
 close to the actual layer sprite image. However, if proactive measures are
 needed to ensure the sprite images exhibit an acceptable resolution, it is
 possible to monitor the bias and gain characteristics of the camera or
 cameras, and subject the images to a correction procedure to compensate
 for any disparity between the camera characteristics. Such corrective
 procedures are well known and so will not be described in detail in this
 specification. Once the images have been corrected, the blending process
 can proceed as described above.
 One unfortunate effect of the blending in Equation (8) is that most forms
 of averaging tend to increase image blur, particularly when there are a
 large number of images. Part of the cause is non-planarity in the scene
 (which will be modeled in the next section), however image noise and
 resampling error also contribute. Therefore, this blurring effect should
 be compensated for (step 504) before proceeding to the residual depth
 estimation module (if employed). One simple method of compensating for
 this effect is to "deghost" the sprite images as described in reference
 [SS98]. Another solution is to explicitly use an image enhancement
 technique such as one of those described in references [IP92, MP94, CB97].
 These latter techniques have the added advantage of estimating a sprite
 with greater resolution than that of the input images.
 1.4 Estimation of Residual Depth
 In general, the scene will not be exactly piecewise planar. To model any
 non-planarity, it can be assumed that the point u.sub.l on the plane
 n.sub.l.sup.T X=0 is displaced slightly. Referring to FIG. 6, it will be
 assumed the point is displaced in the direction of the ray through u.sub.l
 defined by the camera matrix Q.sub.l, and that the distance it is
 displaced is Z.sub.l (u.sub.l), measured in the direction normal to the
 plane (step 600). In this case, the homographic warps used previously are
 not applicable. However, using a similar argument to that in Section 1.2
 and 1.3, it can be shown (see also [KAH94, DTM96]) that:
EQU u.sub.k =H.sup.l.sub.k u.sub.l +Z.sub.l (u.sub.l)t.sub.kl (9)
 where H.sup.l.sub.k =P.sub.k ((n.sub.l.sup.T q.sub.l)I-q.sub.l
 n.sub.l.sup.T)Q.sub.l * is the planar nomography of Section 1.3, t.sub.kl
 =P.sub.k q.sub.l is the epipole, and it is assumed that the plane equation
 vector n.sub.l =(n.sub.x, n.sub.y, n.sub.z, n.sub.d).sup.T has been
 normalized so that n.sup.2.sub.x+ n.sup.2.sub.y+ n.sup.2.sub.z =1.
 Equation (9) can be used to map plane coordinates u.sub.l backwards to
 image coordinates u.sub.k, or to map the image M.sub.kl forwards onto the
 plane (step 602). The result of this warp is denoted by (H.sup.l.sub.k,
 t.sub.kl, Z.sub.l) o M.sub.kl, or W.sup.l.sub.k o M.sub.kl for more
 concise notation.
 To compute the residual depth map Z.sub.l, an optimization of the same (or
 a similar) consistency metric as that used in Section 1.2 to estimate the
 plane equation could be employed (steps 604 and 606). Doing so is
 essentially solving a simpler stereo problem. In fact, almost any stereo
 algorithm could be used to compute Z.sub.l. The one property the algorithm
 should have is that it favors small disparities.
 1.5 Re-Estimation of Layer Sprite Images
 Once the residual depth offsets have been estimated, the layer sprite
 images should be re-estimated using:
 ##EQU2##
 rather than Equation (8). This re-estimation of the layer sprite images
 using the just computed residual depth offsets provides a more accurate
 representation of the per-pixel color values of each sprite image. When
 the images I.sub.k are warped onto the estimated plane of a layer, points
 actually on the plane will register exactly. However, corresponding image
 points in the images that are offset from the estimated plane will warp
 onto the plane in different locations, thereby reducing the accuracy of
 the blended layer sprite image. Thus, re-estimating the layer sprite
 images by first compensating for the parallax effect between points offset
 from the estimated plane and so bringing these points into alignment in
 the warped images, will increase the resolution of the layer sprite image
 estimates. Equation (10) accomplishes this alignment and blending process.
 1.6 Re-Assiqnment of Pixels to Layers
 The initial assignments of pixels to a particular layer were made using one
 of the methods previously described in section 1.1. However, the initial
 assignment process represented a somewhat crude estimate of what layer
 each pixel should belong and was done simply to provide a starting point
 for estimating the plane equations, sprite images and residual depths
 (i.e., n.sub.l, L.sub.l, and Z.sub.l, respectively). The Boolean masks
 B.sub.kl contain this initial assignment information and allowed the
 computation of the masked image M.sub.kl using M.sub.kl
 =B.sub.kl.multidot.I.sub.k. It will now be described how to re-estimate
 the pixel assignments from the estimates of n.sub.l, L.sub.l and Z.sub.l
 to obtain a more accurate grouping of the pixels to the various layers.
 One possible method would be to update the pixel assignments by comparing
 the warped images W.sup.l.sub.k o M.sub.kl to the layer sprite images
 L.sub.l. However, if these images were compared, it would not be possible
 to deduce anything about the pixel assignments outside of the current
 estimates of the masked regions. Thus, it is preferable to make the
 Boolean mask B.sub.kl "grow", thereby allowing for pixel assignments
 outside the current estimated masked regions to be considered. Referring
 to FIG. 7, the Boolean mask B.sub.kl is allowed to "grow" by comparing
 W.sup.l.sub.k o I.sub.k with:
 ##EQU3##
 where Z.sub.l is enlarged so that it declines to zero outside the masked
 region (step 700). The danger with working with warped versions of the
 complete input images W.sup.l.sub.k o I.sub.k is that unoccluded pixels
 will be blended together with occluded pixels and result in poor estimates
 of the "full" layer sprites Full L.sub.l. One possible solution to this is
 a use of a robust blending operator, such as the median. Another solution
 might be to weight masked pixels more than unmasked pixels during the
 blend. A final possibility would be to grow the current estimate of
 B.sub.kl by a small prescribed number of pixels (e.g., 5 pixels), and then
 recompute W.sup.l.sub.k o M.sub.kl. Such a strategy would still allow the
 layer masks to grow slightly on each iteration, but, it is believed, limit
 the influence of occluded pixels in the blend.
 Given the full layer sprites Full L.sub.l, a preferred approach to pixel
 re-assignment is as follows. First compute a measure P.sub.kl (u.sub.l) of
 the likelihood that the pixel W.sup.l.sub.k o I.sub.k (u.sub.1) is the
 warped image of the pixel u.sub.l in the full layer sprite Full L.sub.l
 (step 702). Next, P.sub.kl is warped back into the coordinate system of
 the input image I.sub.k (step 704) to yield:
EQU Warped P.sub.kl =(W.sup.l.sub.k).sup.-1 o P.sub.k (12)
 This warping tends to blur P.sub.kl, but the blurring is believed to be
 acceptable since it is desired to smooth the pixel assignment anyway. In
 fact, it may be advantageous to smooth P.sub.kl even more by, for example,
 using an isotropic smoother such as a Gaussian. Other, more directed,
 smoothing methods could include performing a color segmentation of each
 input image and only smoothing within each segment as described in
 [ASB94]. Alternatively, P.sub.kl might be smoothed less in the direction
 of the intensity gradient since strong gradients often coincide with depth
 discontinuities and hence layer boundaries. These latter two methods may
 produce a better result as it is believed they would tend to smooth
 P.sub.kl where it is most needed, rather than in a general fashion as
 would be the case with the aforementioned Gaussian smoothing procedure.
 The pixel re-assignment can now be computed by choosing the best possible
 layer for each pixel (step 706):
 ##EQU4##
 There are a number of possible ways of defining P.sub.kl. Perhaps the
 simplest is the residual intensity difference [SA96]:
EQU P.sub.k1 =.parallel.W.sup.l.sub.k o I.sub.k -Full L.sub.l.parallel. (14)
 Another possibility is the residual normal flow magnitude:
 ##EQU5##
 Locally estimated variants of the residual normal flow have been described
 in [IP92, IRP92, IAH95]. A third possibility would be to compute the
 optical flow between W.sup.l.sub.k o I.sub.k and Full L.sub.l and then use
 the magnitude of the flow for P.sub.kl.
 Once the pixel re-assignment process is complete and updated Boolean masks
 B.sub.kl have been computed via equation (13), the estimation steps
 described in sections 1.2 through 1.6 can be iteratively repeated until
 the change from one iteration to the next fall below a prescribed
 threshold. This prescribed threshold can relate to one or any combination
 of the layer parameters n.sub.l, L.sub.l, Z.sub.l and B.sub.kl. For
 example, if the change in the Boolean mask image B.sub.kl from the
 immediately preceding iteration falls below a certain percentage change
 threshold, it might be assumed that the then current estimates for the
 other parameters (i.e., n.sub.l, L.sub.l, and Z.sub.l) are as accurate as
 possible for the initial estimation phase of the stereo reconstruction
 process. Alternately, one of the other layer parameters, or a combination
 of any or all of the parameters could be put to a similar threshold test.
 In these alternate methods, if the change between iterations in some
 number or all the tested parameters falls below their associated
 prescribed thresholds, then the iterative process would be stopped and the
 then current layer parameters accepted as the final estimate for the
 initial phase of the process. The actual threshold or thresholds employed
 will depend on the application and the degree of accuracy required. The
 number and value of the threshold(s) can be readily determined given the
 resolution desired for the particular application in which the present
 invention is being employed, and so will not be elaborated upon herein.
 2. Layer Refinement by Re-Synthesis
 In this section, it will be shown how the estimates of the layer sprites
 images, as well as the layer plane equations and residual depths (if
 employed), can be further refined, assuming that the opacities
 .alpha..sub.l associated with the sprite images are now represented by
 real values rather than the Boolean values used previously. The process
 begins by formulating a generative model of the image formation process.
 Afterwards, a proposed measure of how well the layers re-synthesize the
 input images is described, and it is shown how this measure can be
 minimized to refine the complete layer sprite estimates.
 2.1 The Image Formation Process
 Referring to FIG. 8, the first phase of the layer refinement process is to
 formulate a forward model of the image formation process using image
 compositing operations such as described in [Bli94], i.e. by painting the
 sprites one over another in a back-to-front order. The basic operator used
 to overlay the sprites is the over operator:
EQU F.circle-w/dot.B.ident.F+(1-.alpha..sub.F)B.sub.1 (16)
 where F and B are the foreground and background sprites, and .alpha..sub.F
 is the opacity of the foreground. Given this, the generative model
 consists of the following steps:
 a. Using the known camera matrices, plane equations, and residual depths,
 warp each layer backwards onto the coordinate frame of image I.sub.k using
 the inverse of the operator in Section 1.4, (step 800).
 This yields the un-warped sprite image:
EQU U.sub.kl =(W.sup.l.sub.k).sup.-1 o L.sub.l (17)
 Note that the real valued opacities are warped along with the color values.
 b. Composite the un-warped sprite images in back-to-front order:
 ##EQU6##
 to obtain the synthesized image S.sub.k. If the stereo reconstruction
 problem has been solved, S.sub.k should match the input I.sub.k.
 This last step can be re-written as three simpler steps:
 (i) Compute a visibility V.sub.kl for each un-warped sprite image (step
 802) [SG97]:
 ##EQU7##
 where .alpha..sub.kl is the alpha channel of U.sub.kl, and V.sub.kl =1.
 (ii) Compute the masked layer images, M.sub.kl =V.sub.kl U.sub.kl (step
 804).
 (iii) Sum up the masked images, S.sub.k =.SIGMA..sup.L.sub.l=1 M.sub.kl
 (step 806).
 In these three sub-steps, the visibility map makes the contribution of each
 sprite pixel to each image explicit.
 2.2 Minimization of Re-Synthesis Error
 As mentioned previously, if the reconstructed layer estimates are accurate,
 the synthesized image S.sub.k should be similar to the input image
 I.sub.k. If not, the second phase of the layer refinement process is
 employed. Specifically, the prediction error:
 ##EQU8##
 can be used to refine the layer estimates (step 808) using a simple
 technique such as gradient descent. However, rather than trying to
 optimize over all of the parameters (L.sub.l, n.sub.l and Z.sub.l)
 simultaneously, it is preferable to only adjust the sprite colors and
 opacities L.sub.l, and then re-run the previous estimation steps described
 in sections 1.2 through 1.5 to refine n.sub.l and Z.sub.l (step 810).
 In solving Equation (20), the derivatives of the cost function C with
 respect to the colors and opacities in L.sub.l (u.sub.l) can be computed
 using the chain rule [SG97]. In more detail, the visibility map V.sub.kl
 mediates the interaction between the un-warped sprite U.sub.kl, and the
 synthesized image S.sub.k1 and is itself a function of the opacities in
 the un-warped sprites U.sub.kl. For a fixed warping function
 W.sup.l.sub.k, the pixels in U.sub.kl are linear combinations of the
 pixels in sprite L.sub.l. This dependence can either be exploited directly
 using the chain rule to propagate gradients, or alternatively the
 derivatives of C with respect to U.sub.kl can be warped back into the
 reference frame of L.sub.l [SG97].
 While the invention has been described in detail by specific reference to
 preferred embodiments, it is understood that variations and modifications
 thereof may be made without departing from the true spirit and scope of
 the invention.
 References
 [ASB94] S. Ayer, P. Schroeter, and J. Bigun. Segmentation of moving objects
 by robust parameter estimation over multiple frames. In 3rd ECCV, pages
 316-327, 1994.
 [BAHH92] J. R. Bergen, P Anandan, K. J. Hanna, and R. Hingorani.
 Hierarchical model-based motion estimation. In 2nd ECCV, pages 237-252,
 1992.
 [BJ96] M. J. Black and A. D. Jepson. Estimating optical flow in segmented
 images using variable-order parametric models with local deformations.
 PAMI, 18(10):972-986, 1996.
 [Bli94] J. F. Blinn. Jim Blinn's corner: Compositing, part 1: Theory. IEEE
 Computer Graphics and Applications, 14(5):83-87, September 1994.
 [CB97] M. C. Chiang and T. E. Boult. Local blur estimation and
 super-resolution. In CVPR '97, pages 821-826, 1997.
 [DP95] T. Darrell and A. P. Pentland. Cooperative robust estimation using
 layers of support. PAMI, 17(5):474-487, 1995.
 [DTM96] P. E. Debevec, C. J. Taylor, and J. Malik. Modeling and rendering
 architecture from photographs: A hybrid geometry and image-based approach.
 In SIGGRAPH '96, pages 11-20, 1996.
 [IAH95] M. Irani, P. Anandan, and S. Hsu. Mosaic based representations of
 video sequences and their applications. In 5th ICCV, pages 605-611, 1995.
 [IP92] M. Irani and S. Peleg. Image sequence enhancement using multiple
 motions analysis. In CVPR '92, pages 216-221, 1992.
 [IRP92] M. Irani, B. Rousso, and S. Peleg. Detecting and tracking multiple
 moving objects using temporal integration. In 2nd ECCV, pages 282-287,
 1992.
 [KAH94] R. Kumar, P. Anandan, and K. Hanna. Direct recovery of shape from
 multiple views: A parallax based approach. In 12th ICPR, pages 685-688,
 1994.
 [MP94] S. Mann and R. W. Picard. Virtual bellows: Constructing high quality
 stills from video. In 1st ICIP, pages 363-367, 1994.
 [SA96] H. S. Sawhney and S. Ayer. Compact representations for videos
 through dominant and multiple motion estimation. PAMI, 18(8):814-830,
 1996.
 [SG97] R. Szeliski and P. Golland. Stereo matching with transparency and
 matting. In ICCV '98, pages 517-524, 1998.
 [SS97] R. Szeliski and H. -Y. Shum. Creating full view panoramic image
 mosaics and texture-mapped models. In SIGGRAPH '97, pages 251-258, 1997
 [SS98] H. -Y. Shum and R. Szeliski. Construction and refinement of
 panoramic mosaics with global & local alignment. In ICCV '98, pages
 953-958, 1998.
 [WA93] J. Y. A. Wang and E. H. Adelson. Layered representation for motion
 analysis. In CVPR '93, pages 361-366, 1993.