Method for recovering 3D scene structure and camera motion from points, lines and/or directly from the image intensities

An algorithm for recovering structure and motion from points, lines and/or image intensities. The algorithm combines feature based reconstruction and direct methods. The present invention is directed to a method for recovering 3D scene structure and camera motion from image data obtained from a multi-image sequence, wherein a reference image of the sequence is taken by a camera at a reference perspective and one or more successive images of the sequence are taken at one or more successive different perspectives by translating and/or rotating the camera. The method comprising the steps of (a) determining image data shifts for each successive image with respect to the reference image; the shifts being derived from the camera translation and/or rotation from the reference perspective to the successive different perspectives; (b) constructing a shift data matrix that incorporates the image data shifts for each image; (c) calculating two rank-3 factor matrices from the shift data matrix using SVD, one rank-3 factor matrix corresponding the 3D structure and the other rank-3 factor matrix corresponding the camera motion; (d) recovering the 3D structure from the 3D structure matrix using SVD by solving a linear equation; and (e) recovering the camera motion from the camera motion matrix using the recovered 3D structure.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention is directed to a method for recovering 3D structure and camera motion, and more particularly to a linear algorithm for recovering the structure and motion data from points, lines, and/or directly from the image intensities.

2. Prior Art

The science of rendering a 3D model from information derived from a 2D image predates computer graphics, having its roots in the fields of photogrammetry and computer vision.

Photogrammetry is based on the basic idea that when a picture is taken, the 3D world is projected in perspective onto a flat 2D image plane. As a result, a feature in the 2D image seen at a particular point actually lies along a particular ray beginning at the camera and extending out to infinity. By viewing the same feature in two different photographs the actual location can be resolved by constraining the feature to lie on the intersection of two rays. This process is known as triangulation. Using this process, any point seen in at least two images can be located in 3D. It is also possible to solve for unknown camera positions as well with a sufficient number of points. The techniques of photgrammetry and triangulation were used in such applications as creating topographic maps from aerial images. However the photogrammetry process is time intensive and inefficient.

Computer vision techniques include recovering 3D scene structure from stereo images, where correspondence between the two images is established automatically from two images via an iterative algorithm, which searches for matches between points in order to reconstruct a 3D scene. It is also possible to solve for the camera position and motion using 3D scene structure from stereo images.

Current computer techniques are focused on motion-based reconstruction and are a natural application of computer technology to the problem of inferring 3D structure (geometry) from 2D images. This is known as Structure-from-Motion. Structure from motion (SFM), the problem of reconstructing an unknown 3D scene from multiple 2D images, is one of the most studied problems in computer vision.

Most SFM algorithms that are currently known reconstruct the scene from previously computed feature correspondences, usually tracked points. Other algorithms are direct methods that reconstruct from the images intensities without a separate stage of correspondence computation. Previous direct methods were limited to a small number of images, required strong assumptions about the scene, usually planarity or employed iterative optimization and required a starting estimate.

These approaches have complementary advantages and disadvantages. Usually some fraction of the image data is of such low quality that it cannot be used to determine correspondence. Feature-based method address this problem by pre-selecting a few distinctive point or line features that are relatively easy to track, while direct methods attempt to compensate for the low quality of some of the data by exploiting the redundancy of the total data. Feature-based methods have the advantage that their input data is relatively reliable, but they neglect most of the available image information and only give sparse reconstructions of the 3D scene. Direct methods have the potential to give dense and accurate 3D reconstructions, due to their input data's redundancy, but they can be unduly affected by large errors in a fraction of the data.

A method based on tracked lines is described in “A Linear Algorithm for Point and Line Based Structure from Motion”, M. Spetsakis, CVGIP 56:2 230-241, 1992, where the original linear algorithm for 13 lines in 3 images was presented. An optimization approach is disclosed in C. J. Taylor, D. Kriegmann, “Structure and Motion from Line Segements in Multiple Images,” PAMI 17:11 1021-1032, 1995. Additionally, in “A unified factorization algorithm for points, line segments and planes with uncertainty models” K. Morris and I. Kanade, ICCV 696-702, 1998, describes work on lines in an affine framework. A projective method for lines and points is described in “Factorization methods for projective structure and motion”, B. Triggs, CVPR 845-851, 1996, which involves computing the projective depths from a small number of frames. “In Defense of the Eight-Point Algorithm: PAMI 19, 580-593, 1995, Hartley presented a full perspective approach that reconstructs from points and lines tracked over three images.

The approach described in M. Irani, “Multi-Frame Optical Flow Estimation using Subspace Constraints,” ICCV 626-633, 1999 reconstructs directly from the image intensities. The essential step of Irani for recovering correspondence is a multi-frame generalization of the optical-flow approach described in B. Lucas and T. Kanade, “An Iterative Image Registration Technique with an Application to Stereo Vision”, IJCAI 674-679, 1981, which relies on a smoothness constraint and not on the rigidity constraint. Irani uses the factorization of D simply to fill out the entries of D that could not be computed initially.

Irani writes the brightness constancy equation (7) in matrix form as Δ=−DI, where D tabulates the shifts diand I contains the intensity gradients ∇I(pn). Irani notes that D has rank 6 (for a camera with known calibration), which implies that Δ must have rank 6. To reduce the effects of noise, Irani projects the observed Δ onto one of rank 6. Irani then applies a multi-image form of the Lucas-Kanade approach to recovering optical flow which yields a matrix equation DI2=−Δ2, where the entries of I2are squared intensity gradients IaIbsummed over the “smoothing” windows, and the entries of Δ2have the form IaΔI. Due to the added Lucas-Kanade smoothing constraint, the shifts D or dnican be computed as D=−Δ2[I2]+denotes the pseudo-inverse, except in smoothing windows where the image intensity is constant in at least one direction. Using the rank constraint on D, Irani determines additional entries of D for the windows where the intensity is constant in one direction.

SUMMARY OF THE INVENTION

The present invention is directed to a method for recovering 3D scene structure and camera motion from image data obtained from a multi-image sequence, wherein a reference image of the sequence is taken by a camera at a reference perspective and one or more successive images of the sequence are taken at one or more successive different perspectives by translating and/or rotating the camera. The method comprising the steps of:(a) determining image data shifts for each successive image with respect to the reference image; the shifts being derived from the camera translation and/or rotation from the reference perspective to the successive different perspectives;(b) constructing a shift data matrix that incorporates the image data shifts for each image;(c) calculating two rank-3 factor matrices from the shift data matrix using SVD, one rank-3 factor matrix corresponding the 3D structure and the other rank-3 factor matrix corresponding the camera motion;(d) recovering the 3D structure from the 3D structure matrix using SVD by solving a linear equation; and(e) recovering the camera motion from the camera motion matrix using the recovered 3D structure.

The method of the invention is a general motion algorithm wherein the camera positions for each successive perspective do not lie on a single plane. The method can reconstruct the image from points, lines or intensities.

The present invention is an essentially linear algorithm that combines feature-based reconstruction and direct methods. It can use feature correspondences, if these are available, and/or the image intensities directly. Unlike optical-flow approaches such as that described in “Determining Optical Flow”, B. K. P Horn and B. G. Schunck, AI 17, 185-203, 1981, the algorithm exploits the rigidity constraint and needs no smoothing constraint in principle assuming the motion is small and the brightness constancy equation is valid. However, in practice, it is sometimes necessary to impose smoothness in the algorithm.

The method of the present invention can reconstruct from both point and line features. This is an important aspect of the present invention, since straight lines are abundant, especially in human-made environments, and it is often possible to track both feature types. Lines can be localized very accurately due to the possibility of finding many sample points on a single line, and they do not suffer from some of the problems of point tracking such as spurious T-junctions and the aperture effect. The method described here is the first that can reconstruct from all three types of input data in any combination or separately over a sequence of arbitrarily many images under full perspective.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

Definitions

For purposes of this disclosure, the invention will be described in accordance with the following terms defined herein. The method of the present invention assumes that the calibration is known and takes the focal length to be 1 (wlog). Modifications to the method for uncalibrated sequences will be discussed later.

Assuming a sequence of Niimages. Choose the zeroth image as the reference image. Let Ri, Tidenote the rotation and translation between this image and image i. Parameterize small rotations by the rotational velocity ω≡(ωxi,ωyi,ωzi)T. Let [Ti]2≡(TxiTyi)T.

Let a 3D point P transform as P′=R(P−T). Assume Nppoints and NLlines tracked over the sequence. For clarity, we use a different notation for the tracked points than for the pixel positions in the intensity images. Let qmi≡(qx,qy)midenote the m-th tracked point and Ali≡(α,β,γ)lTidenote the i-th line in the i-th image. Let {overscore (q)}≡(q;1)T. Let R*q denote the image point obtained from q after a rotation: R*q≡[(R{overscore (q)})x,(R{overscore (q)})y]T/(R{overscore (q)})z. Define the three point rotational flows r(1)(x,y), r(2)(x,y), r(3)(x,y), by[r(1),r(1),r(1)]≡[(-xy-(1+y2)),(1+x2xy),(-yx)].
{circumflex over (q)}midenotes the 3D unit vector {overscore (q)}/|{overscore (q)}|=(q;1)T/|(q;1)|. Similarly, Â=A/|A|. Let smi≡qmi=qm0denote the image displacement for the m-th tracked point and Qm=(Qx, Qy, Qz)mTdenote the 3D scene point corresponding to the tracked point qm. Then parameterize a 3D line L by two places that it lies in. The first plane, described by A, passes through the center of the projection and the imaged line in the reference image. The normal to the second plane, B, is defined by requiring B·A=0 and B·Q=−1 for any point Q on L.

Let pn≡(xn·yn)Tbe the image coordinates of the n-the pixel position. Order the pixels from 1 to Nx, the total number of pixels. Let Iidenote the i-the intensity image and let Ini=Ii(pn) denote the image intensity at the n-the pixel position in Il. Let Pndenote the 3D point imaged at pnin the reference image, with Pn≡(Xn, Xn, Xn)Tin the coordinate system of I0. Let dnidenote the shift in image position from Ioto Iiof the 3D point Pn.

Suppose Vαis a set of quantities indexed by the integer α. We use the notation {V} to denote the vector with elements given by the Vα.

Algorithm Description-Preliminary Processing

The method of the present invention requires that the translational motion not be to large, (e.g., with |T/Zmin|≦1/3 and that the camera positions do not lie in a plane.

Given tracked points and lines, we approximately recover the rotations between the reference image and each following image by minimizing:∑m⁢(q^ni-Ri⁢qmo)2+μ⁢∑l⁢(A^li-Ri⁢A^l)2(1)
with respect to the rotations Ri, where one should adjust the constant according to the relative accuracy of the point and line measurements. Minimizing equation (1) gives accurate rotation as long as the translations are not too large. After recovering the rotations, compensation is made for them and henceforth they can be treated as small.
Points

Point displacements are calculated in accordance with J. Oliensis, “A Multi-Frame Structure from Motion Algorithm under Perspective Projection”, IJCV 34:2/3 163-192, 1999, incorporated herein by reference, the point displacements aresmi=Qz,m-1⁡(qm⁢Tzi-[Ti]2)1-Qz,m-1⁢Tzi+f⁡(Ri,qmi)
where the rotational flow ƒ(Ri,qmi)≡qmi−(Ri)−1*qmi. Under the assumptions, we can approximate smi≈Qz,m−1(qmTzi−[Ti]2)+ωxir(1)(qm)+ωyxir(2)(qm)+ωzir(3)(qm) where the last three terms give the first-order rotational flow. Correspondingly, define the three length-2Nptranslational flow vectorsΦx≡-[(Qz-1}{0}]⁢⁢Φy≡-[{0}{Qz-1}]⁢⁢Φz≡-[{qx⁢Qz-1}{qy⁢Qz-1}]
and the length −2Nprotational flowsΨx=[{rx(1)⁡(q)}{ry(1)⁡(q)}]⁢⁢Ψy≡[{rx(2)⁡(q)}{ry(2)⁡(q)}]⁢⁢Ψz≡[{rx(3)⁡(q)}{ry(3)⁡(q)}]

Then collect the shifts, smismiinto a Ni×2Npmatrix S, where each row corresponds to a different image i and equals [{sxi}T{syi}T]. Then
S≈{Tx}ΦxT+{Tx}ΦxT+{Tx}ΦxT+{ωx}ΨxT+{ωy}ΨyT+{ωz}ΨzT.  (2)
Lines

The line measurement model is described as follows. One can reasonably assume that the noise of a measured line is proportional to the integrated image distance between the measured and true positions of the line. Let ΘFOVdenote the field of view (FOV) in radians. Typically, ΘFOV<1. If the measured line differs from the true one by a z-rotation, this gives a noise of O((ΘFOV/2)2). A rotation around an axis in the x-y plane gives a noise of O(ΘFOV), which is typically larger. (These estimates reflect the fact that the displacement region bounded by the true and measured lines look like two triangles in the first case and a quadrilateral in the second. For the representation described herein, this implies that line measurements determine Azmore accurately than Ax, Ayby a factor roughly of 1/ΘFOV.

The Image flow for lines can be described as follows. Let A≡(A;0) and B≡(B;−1) be the homogeneous length-4 vectors corresponding to the plane normals A and B for the line L. Let Q≡(Q;1) be the homogeneous vector corresponding to a point on L. Then A·B=A·Q=B·Q=0. After rotation R and translation T, we determine the transformed A′B′ from the requirement that A′·Q′=B′·Q′=0. We can satisfy this requirement usingA*=[R⁢⁢AT·A],B*=[R⁢⁢BT·B+B4].
But A* doesn't necessarily have A4′=0, as the A′ for the new image of the line must. So we set A′≡A*−A4*B*/B4*, which implies that the new image line is given byA′≡R⁡(A-B⁢T·AT·B+B4).(3)
Now consider the line flow δAli=Ali−Al0. Define δAR=Al−R−1Aiandδ⁢⁢AT≡R-1⁢Ai-A=(T·A)1-T·B⁢B⁢⁢so⁢⁢δ⁢⁢A=δ⁢⁢AT+δ⁢⁢AR.
For small rotations and translations, with |T·B|<<1,
δAi≈(T·A)B+ωΔA.(4)B·A=0 implies that δA·A≈0. To eliminate the scale ambiguity in defining each A, we take |A10|=1 in the reference image and require that the line flow satisfies 0=δAli·Alo≡(Ali−Al0)·Aloexactly. We normalize all A0to the same magnitude in the reference image to reflect the fact that the measurement errors should be similar for all lines. The requirement that 0=δAli·Alo“fixes the gauge” in a way that is consistent with our line representation and small-motion assumption.

After “gauge fixing” each δAliincorporates just two degrees of freedom. We represent δAliby its projection along two directions, Al×({circumflex over (z)}×Al) and {circumflex over (z)}×Al, which we refer to respectively as the upper and lower directions. Let the unit 3-vector PUland PLlproject onto these two directions. For typical ΘFOV<1,|Al·{circumflex over (z)}|<<1 and the upper component of Alroughly equals Al,z. Thus, image measurements determine the upper component more accurately than the lower, by roughly 1/ΘFOV. In analogy to the point definitions, define the line translation flowsΦL⁢⁢x≡[{Ax⁢BU}{Ax⁢BL}]⁢⁢ΦL⁢⁢y≡[{Ay⁢BU}{Ay⁢BL}]⁢⁢Φz≡[{Az⁢BU}{Az⁢BL}]
where these are length-2NLvectors. BU≡B·PUand BL≡B·PLare the upper and lower components of B. Similarly, define the rotational flowsΨL⁢⁢x≡[{PU·(x×A)}{PL·(x×A)}]⁢⁢ΨL⁢⁢y≡[{PU·(y×A)}{PL·(y×A)}]⁢⁢Ψz≡[{PU·(z×A)}{PL·(z×A)}](5)

Let Λ be the N1×2NLmatrix where each row corresponds to a different image i and equals [{PU·δAi}T{PU·δAi}T. Then
Λ≈{Tx}ΦLxT+{Ty}ΦLyT+{Tz}ΦLzT+{ωyz}ΨLxT+{ωyzT}ΨLyT+{ωzz}ΨLzT(6)
Intensities

One can apply the same arguments to the dni, the image shifts corresponding to the 3D point imaged at the pixel position pn, as to the smifor the tracked feature points. Thus dni≈Zn−1(pnTzi−[Ti]2)+ωxir(1)(pn)+ωxir(2)(pn)+ωxir(3)(pn). Let ∇In≡∇I(pn)≡(Ix,Iy)nTrepresent the gradient of the image intensity for I0, with some appropriate definition for a discrete grid of pixels. Let ΔInidefine the change in intensity with respect to the reference image. Its simplest definition is ΔIni=Ini−In0. The brightness constancy equation is
ΔIni+∇In·dni=0.  (7)

The reconstruction of the 3D scene is described as follows. Define the total rotational flow vectors for points, lines and image intensities by {overscore (Ψ)}xT≡[ΨxTωLΨlxTωlΨlxT], and similarly for the y and z components. Here ωLand ω1are constant weights that the user should set according to the relative noisiness of the point, line, and intensity data. The {overscore (Ψ)}αhave length 2Np+2NL+NX≡Ntot. One can verify from the definitions that the {overscore (ω)} are computable from measured quantities and can be taken as known. Using Householder matrices, one can compute a (Ntot−3)×Ntotmatrix H annihilating the three {overscore (ω)}x,y,z. Computing and multiplying by H cost O(Ntot) computation. Define the N1×Ntotmatrix {overscore (D)}≡[S ωLA ωLΔ]. Let C be a constant (N1−1)×(N1−1) matrix with Cii′≡δii′+1, we include C to counter the bias due to singling out the reference image for special treatment. DefineD_≡C-12⁢D_⁢⁢HT.
Define the total translational flow vectors byΦ_x≡[ΦxωL⁢ΦLxω1⁢Φ1⁢x],
and similarly for y and z. From above, equations (2), (6) and (8) imply thatD_CH≈C-12⁢{Tx}⁢Φ_T⁢HT+C-12⁢{Ty}⁢Φ_yT⁢HT+C-12⁢{Tz}⁢Φ_zT⁢HT.(9){overscore (D)}CHis approximately rank 3.

Our basic algorithm of the method of the present invention is,1) Define H and compute {overscore (D)}CH. Using the singular value decomposition (SVD), compute the best rank-3 factorization of {overscore (D)}CH≈M(3)S(3)Twhere M(3), S(3)are rank 3 matrices corresponding respectively to the motion and structure.2) S(3)satisfies,
[{overscore (Φ)}x{overscore (Φ)}y{overscore (Φ)}z]=HTS(3)U+[{overscore (Ψ)}x{overscore (Ψ)}y{overscore (Ψ)}z]Ω  (10)
where U and Ω are unknown 3×3 matrices. We eliminate the unknowns Qz, Bzand Z−1from the {overscore (Φ)}αin this system of equations to get 3Ntotlinear constraints on the 18 unknowns in U and Ω. We solve these constraints with O(Ntot) computations using the SVD.3) Given U and Ω, we recover the structure unknowns Qz, Bzand Z−1from [{overscore (Φ)}x{overscore (Φ)}x{overscore (Ψ)}x]=HTS(3)U+[{overscore (Ψ)}x{overscore (Ψ)}y{overscore (Ψ)}z]Ω where Q is the 3d coordinate for an image pixel in the reference image, B is the shortest vector from the camera center to a 3D line and Z is the depth from the camera to a 3D scene alone the cameras optical axis.4) Given U, we use S(3)U≈[{overscore (Φ)}x{overscore (Φ)}y{overscore (Φ)}z] andD_CH≈C-12⁢{Tx}⁢Φ_T⁢HT+C-12⁢{Ty}⁢Φ_yT⁢HT+C-12⁢{Tz}⁢Φ_zT⁢HT
to recover the translations.5) We recover the rotations ωxi, ωyi, ωzifromωxi⁢Ψ_xn+ωyxi⁢Ψ_yn+ωzx⁢Ψ_zn=C-12⁢D_ni-(C-12⁡({Tx}⁢Φ_T+{Ty}⁢Φ_yT+{Tz}⁢Φ_zT))ni

The description above omits some steps, which are important for bias correction. 1) The upper line components in {overscore (D)} are weighted by a factor proportional to 1/ΘFOV, to account for the greater accuracy in the measurement of these components. 2) To correct for bias due to the fact that the FOV is typically small, in Steps2-4we weight Tziby a factor proportional to the FOV and weight {overscore (Φ)}zby the inverse of this, while leaving the x and y Components untouched. 3) The steps of the algorithm can be iterated to give better results and correct for the small motion assumptions. This involves multiplying the original feature point shifts by 1−Z−1Tzand the line shifts by 1−B·T. The algorithm is guaranteed to converge to the correct reconstruction if the motion and noise are small and the camera positions do not lie on a plane.

Of course, the results for the intensity-based part of our algorithm depend crucially on the technique used for computing derivatives. To make this more robust, one should iteratively reduce the size of the displacements dniby warping and recomputing the reconstruction in a coarse-to-fine mode. One implementation simply computes the image derivatives at a single scale, using consistency between the spatial derivatives computed for different images to determine whether the assumption of small image motion holds for this current scale. Only those pixel sites where the assumption does hold are used for the motion computation.

We have sometimes found it useful to preprocess the image intensities prior to running the algorithm. We first compute a Laplacian image and then transform the intensities by a stigma function to enhance edgy regions and suppress textureless ones. This gives the intensity image a form that is intermediate between the unprocessed images and a selected set of tracked features.

Bas-Relief for Lines

Due to the bas-relief ambiguity, it is difficult to recover the constant component of the vector {Z−1} for point features, though typically one can recover all other components accurately. A similar result can be derived for lines. The rotational flow for a line is ω×A. For typical ΘFOV<<1,|Az|<<|Ax,y|. Thus {circumflex over (x)} and ŷ rotations give flows roughly proportional to ŷ×A˜−Ax{circumflex over (z)} and ({circumflex over (x)}×A)˜−Ay{circumflex over (z)}. From equation (4) above, the effects of a z-translation Tzare suppressed by |Az|, and the translational flows due to the constant component of {Bz} constant component from rotational effects, which makes this component difficult to recover accurately.

One can get a more quantitative analysis of this effect as follows. Assuming that the input data consists of lines only, [{overscore (Φ)}x{overscore (Φ)}x{overscore (Φ)}x]=HTS(3)U+[{overscore (Ψ)}x{overscore (Ψ)}y{overscore (Ψ)}z]Ω implies that {overscore (P)}sH[{overscore (Φ)}Lx, {overscore (Φ)}Ly, {overscore (Φ)}Lz]=0 where {overscore (P)}s,is a projection matrix that annihilates S(3). Those {B} components that produce small overlaps∑a⁢Φ_LaT⁢HT⁢Ps⁢H⁢⁢Φ_la,
will also produce small overlaps∑a⁢Φ_LaT⁢HT⁢Ps⁢H⁢⁢Φ_La.
It follows from standard perturbation theory that noise will mix these components into the true solution for the B, where the amount of contamination is inversely proportional to the size of the eigenvalue. Define B≡[{Bx}; {By}; {Bz}]. We computed for several sequences the eigenvalues of HB, defined such thatBT⁢HB⁢B=∑a⁢Φ_LaT⁢HT⁢H⁢⁢Φ_La,
and found that, as expected, there was one small eigenvalue, whose eigenvector approximately equaled the constant component of {Bz}. Thus one can use the input data to estimate the inaccuracy in recovering this component.
General-Motion Algorithm for Intensities

The General Motion Algorithm for just Intensities is as follows.0. First, assuming that the translations are zero, we recover the rotations and warp all images I1. . . INl−1toward the reference image I0. In this case we then let the image displacements dnirefer to the unrotated images.1. We then compute H andΔCH≡C-12⁢Δ⁢⁢HT,
and using the singular value decomposition (SVD), compute the best rank-3 factorization of −ΔCH{overscore (D)}CH≈M(3)S(3)T, where M(3)and S(3)are rank 3 and correspond respectively to the motion and structure.2a. Further S(3)satisfies, Φ=HTS(3)U+ΨΩ
where U and Ω are unknown 3×3 matrices. We eliminate the unknowns Z−1from the {overscore (Φ)}αin this system of equations to get 3Nplinear constraints on the 18 unknowns in U and Ω. We solve these constraints with O(Np) computations using the SVD. Then given U and Ω, we recover the structure unknowns Z−1from (5).2b. Then given U, we use S(3)U≈Φ and equation (4) to recover the translations.3. The rotations W can then be recovered from-C-12⁢W⁢⁢ΨT=-C-12⁡(T_⁢⁢ΦT+Δ).
Algorithm Analysis

Our neglect of the translations in the basic general-motion algorithm above produces small errors in recovering the rotations when the translations are moderate, with |Ttruei|/Zmin≦1/4 as described in “Rigorous Bounds for Two-Frame Structure from Motion,” J. Oliensis IUW, 1225-1231, 1994. We describe two multi-frame techniques for recovering the rotations. Neglecting the translations in (4) gives-C-12⁢Δ⁢≈-12⁢W⁢⁢ΨT,
which can be solved for W. This technique resembles the projective algorithm of L. Zelnik-Manor and M. Irani, described in “Multi-Frame Alignment of Planes”, CVPR I:151-156, 1999 for reconstructing planar scenes, except that our use of-C-12
reduces the bias due to overemphasizing the noise in the reference image. As usual, one can extend the technique to larger rotations by iteratively warping and re-solving for W.

The technique just described is best for small rotations, since it uses the first-order rotational displacements. Another approach begins with the extension to planar scenes as described in “Structure from Motion from Points, Lines and Intensities, J. Oliensis and M. Werman CVPR 2000 of the multi-frame Sturm-Triggs technique set forth in “A factorization based algorithm for multi-image projective structure and motion”, P. Sturm and B. Triggs, ECCV 709-720, 1996. This multi-frame approach recovers homographies of arbitrary size between the reference image and the subsequent images. One can recover the rotations by choosing the orthogonal, positive-determinant matrices closest to the recovered homographies.

From (4), the columns of S(3)generate approximately the same subspace as that generated by the Φx, Φy, Φz;
S(3)U≈HΦwhere U is a 3×3 matrix. (7) Above gives an overconstrained system of linear equations that one can solve directly for U and Zn−1. However this is an O(Np3) computation. For large images, we use instead the following O(Np) technique.

We eliminate the Zn−1from (8) above and solve directly for U and Ω. Denote the columns of U by U≡[U1U2U3], and similarly for Ω. Let U3′≡s−1U3and Ω3′≡s−1Ω3, where the scale s equals the average distance of the image points from the image center. We include s to reduce the bias of the algorithm. From the definitions of Φx, Φy, Φz, (8) above implies Iyn[HTS(3)U1+ΨΩ1]n≈Ixn[HTS(3)U2+ΨΩ2]n, ≈Iyn[HTS(3)U3′+ΨΩ3′]n−s−1(∇In·pn)[HTS(3)U2+ΨΩ2]n≧Iyn[HTS(3)U3′+ΨΩ3′]nand a similar equation with (xy). Step 2A of our algorithm solves this system of linear equations for Ω and U in the least-squares sense and then recovers the Zn−1from these solutions and (8) above. The computation is O(Np). Note also that step 2a bases its computations on ∇I. It we use the value of ∇I computed from the measured reference image I0, then the estimates of U, Ω, Zn−1multi-frame estimates. To get a better multi-frame estimate, one can first re-compute I0and ∇I based on all the image data and use the result to compute U, Ωvia (9) as follows.

Let ΔCH(3)be the best rank 3 approximation toΔCH≡C-12⁢Δ⁢⁢HT.
Up to unknown rotational flows, ΔCH(3)H gives the intensity shifts due to the translational displacements:C-12⁢Δ≈ΔCH(3)⁢H+W⁢⁢ΨT

We then solve (10) in the least-squares sense for W, which gives improved estimates Δefor Δ and Ie0for I0. These can then be used in (9). However, this computation of Ie0gives O(Nl−1ω2, Nl−2τ2Z−2) errors. If the original noise in I0is less than this, one should use I0directly in step 2a.

We haveC-12⁢T_≈M(3)⁢UT,
where UTis also a 3×3 matrix. The algorithm recovers UTand {overscore (T)} by solving the linear system −ΔCH≈M(3)UTΦTHT, for UT, using the Φ computed previously, and then plugging in the result to recover {overscore (T)}.

General Notes

Experimentally, we find that the intensity pattern varies significantly from image to image in our sequences. We actually apply our approach to a modified version of the original sequence, obtained by:1) Filtering the sequence with a laplacian of gaussian to emphasize edges;2) Applying a nonlinear transformation (a sigma function) to further emphasize the high-frequency features. The new sequences give a continuous representation of the interest features, which can be used to compute derivatives. This procedure emphasizes the discrete, more easily trackable features, but does not eliminate the information from other image regions. One could extend this approach to deal with a variety of interest features, selecting among them to determine the derivatives that are likely to give the best flow. To ensure that the BCE is approximately valid, we also require that the spatial derivatives be consistent over the whole sequence, and that the image intensity be well approximated by a plane at the scale at which we are working. In our current implementation, we apply the algorithm at a single scale.

In another embodiment of the present invention, the algorithm described above can also handle: cameras with changing and unknown focal lengths, where the calibration is otherwise fixed and known and camera calibrations that change arbitrarily from image to image in an unknown fashion (this is know as the projective case).

For varying, unknown focal lengths, one can modify the algorithm along the line of the method described in “A Multi-frame Structure from Motion Algorithm under Perspective Projection” IJCV 34:2/3, 163-192, 199 and Workshop on Visual Scenes 77-84, 1999. The modification is described as follows, in Step 2a, we define H to annihilate ΨF≡{pn·∇I(pn)} in addition to Ψx, Ψy, Ψz. In recovering the Zn−1one must replace Ψ in (8) by [Ψ, ΨF], and Ω becomes a 4×3 matrix. The modified Step 2a can recover the Z−1except for the constant (Zn=1) component, which is intrinsically difficult to recover accurately when the focal lengths are unknown and varying.

For the projective case, the algorithm can be modified along the lines of the method described in “Fast Algorithms for Projective Multi-Frame Structure from Motion”, J. Oliensis and Y. Genc, ICCV 536-543, 1999. In this case, we define H to annihilate eight length-Npvectors, where the n-th entries of these vectors are given by the eight quantities (∇In)α, Pnb(∇In)c, Pnd(∇In)·pn, for a, b, c, dε{x,y}. Step 2a must be modified by replacing Ψ in (8) by a matrix consisting of these eight vectors.

Projective Algorithm

The algorithm described here for calibrated sequences generalizes in a straightforward way to deal with uncalibrated sequences. In this projective case, instead of a preliminary stage of rotation computation, one computes planar homographies between the reference image and each subsequent image. It is worth noting that one can easily and accurately compute these homographies by a multi-frame generalization of the projective-reconstruction method. We briefly describe how this works for the example of tracked points. Assume the scene is planar. Then
λmiqm−i=MiSm(11)
where Miis a 3×3 matrix (a homography), Smis the structure 3-vector giving the position of the m-th point on the plane (in some arbitrary coordinate system), and λmiis the projective depth. The steps of the homography recovery are:1) Take λmi=1 initially.2) For the current estimate of the λmi, collect the λmi{overscore (q)}miinto a single 3N1×Npmatrix Γ. Use the SVD to decompose this into the product of two rank 3 factors: Γ≈M(3)S(3).3) For the current M(3), S(3), compute the λmiminimizing the |Γ−M(3)S(3)|. Return to Step 2.

After convergence, M(3)contains the Miestimates. Our estimate of the homography taking the reference image to image i is M0\Mi. This procedure is guaranteed to converge to the local minimum ofΓ(4)2Γ2,
where Γ(4)is the difference between Γ and its best rank-3 approximation. Also, it gives nearly the optimal estimates of the Miif the true motion is not too large and the scene is truly planar.

The uncalibrated version of the line-based algorithm is described as follows. Let Kibe the standard 3×3 matrix, with zero lower-diagonal entries, giving the calibration parameters for the i-th image. For the uncalibrated case, equation (3) still holds if one replaces R→(K′)−TRK≡MHand T→KT, where K and K′ are the calibration matrices for the first and second image. MHis a general homography and is the inverse transpose of the analogous homography for points. The first order approximation of the new version ofA′=R⁡(A-B⁢T·AT·B+B4)
is δAi≈(T·A)B+ω×A, but with T→KT and with ω×A replaced by δMHA, where MH≡1+δMH.

We define H to annihilate the first-order homography flow due to δMH, instead of the rotational flows as for the calibrated case. Since the first-order approximation of MH−Tis 1−δMHT, one can easily define H to annihilate the first-order flows for both points and lines. RepresentMH≡[FGJT0],
where F is 2×2 and G and J are 2×1. The first-order point and line displacements due to each of the 8 parameters in δMHare given in the following table:

We define H to annihilate the 8 vectors associated with the columns of this table. After the indicated replacements, the uncalibrated algorithm is the same as the calibrated one. If one fixes the point coordinates in some reference image, the remaining freedom under a projective transform is Z−1→ax+by+c+dZ−1where x and y are the image coordinates. This corresponds to scaling and adding an arbitrary plane to the structure. One can derive a similar relation for B. Let P be a 3D point on the line determined by B. Then B·P=−1 implies B·(x,y,1)=−Z−1. A projective transform must preserve B′·P′=−1. From this, and the transform of Z−1given above, it follows that B′−B=(a, b, c) up to a scaling, with the same constants a, b, c as for the point transform. This relation holds with the same constants for any line B.

The fact that B is recoverable only up to an overall shift is also clear from the uncalibrated version of the algorithm. When we use H to eliminate the first-order effects of small homographies, we also eliminate the translational image flows due to the three constant components of the B. Thus, these components cannot be recovered.

Implementation

It will be apparent to those skilled in the art that the methods of the present invention disclosed herein may be embodied and performed completely by software contained in an appropriate storage medium for controlling a computer.

Referring toFIG. 1, which illustrates in block-diagram form a computer hardware system incorporating the invention. As indicated therein, the system includes a video source101, whose output is digitized into a pixel map by a digitizer102. The digitized video frames are then sent in electronic form via a system bus103to a storage device104for access by the main system memory during usage. During usage the operation of the system is controlled by a central-processing unit, (CPU)105which controls the access to the digitized pixel map and the invention. The computer hardware system will include those standard components well-known to those skilled in the art for accessing and displaying data and graphics, such as a monitor,106and graphics board107.

The user interacts with the system by way of a keyboard108and or a mouse109or other position-sensing device such as a track ball, which can be used to select items on the screen or direct functions of the system.

The execution of the key tasks associated with the present invention is directed by instructions stored in the main memory of the system, which is controlled by the CPU. The CPU can access the main memory and perform the steps necessary to carry out the method of the present invention in accordance with instructions stored that govern CPU operation. Specifically, the CPU, in accordance with the input of a user will access the stored digitized video and in accordance with the instructions embodied in the present invention will analyze the selected video images in order to extract the 3D structure and camera motion information from the associated digitized pixel maps.

Referring now toFIG. 2the method of the present invention will be described in relation to the block diagram. A reference image of the sequence is taken by a camera at a reference perspective and one or more successive images of the sequence are taken at one or more successive different perspectives by translating and/or rotating the camera in step201. Me images are then digitized202for analysis of the 3D image content, i.e. points, lines and image intensities. Then image data shifts for each successive image with respect to the reference image are determined203; the shifts being derived from the camera translation and/or rotation from the reference perspective to the successive different perspectives. A shift data matrix that incorporates the image data shifts for each image is then constructed204and two rank-3 factor matrices from the shift data matrix using SVD, one rank-3 factor matrix corresponding the 3D structure and the other rank-3 factor matrix corresponding the camera motion are calculated205. Recovering the 3D structure from the 3D structure matrix using SVD by solving a linear equation206. The camera motion can then be recovered using the recovered 3D structure207.