# Patent

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

## Abstract

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
RELATED APPLICATIONS
The present application is related to U.S. application Ser. No. 09/652,820, titled A Method for Recovering 3D Scene Structure and Camera Motion Directly from Image Intensities, filed on Aug. 31, 2000, by the same inventor as the present application which related application is incorporated herein by reference. This application claim benefit of Provisional 60/210,099 filed Jun. 7, 2000.


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 di and 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 I2 are squared intensity gradients Ia Ib summed over the “smoothing” windows, and the entries of Δ2 have the form Ia ΔI. Due to the added Lucas-Kanade smoothing constraint, the shifts D or dni can 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.



BRIEF DESCRIPTION OF THE DRAWINGS
These and other features, aspects, and advantages of the methods of the present invention will become better understood with regard to the following description, appended claims, and accompanying drawings where:
FIG. 1 schematically illustrates a hardware implementation of the present invention.
FIG. 2 is a block diagram that illustrates the method of the present invention.



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 Ni images. Choose the zeroth image as the reference image. Let Ri, Ti denote 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 Np points and NL lines 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)mi denote the m-th tracked point and Ali≡(α,β,γ)lTi denote 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
                      +
                      
                        y
                        2
                      
                    
                    )
                  
                
              
            
          
          )
        
        ,
        
          (
          
            
              
                
                  1
                  +
                  
                    x
                    2
                  
                
              
            
            
              
                xy
              
            
          
          )
        
        ,
        
          (
          
            
              
                
                  -
                  y
                
              
            
            
              
                x
              
            
          
          )
        
      
      ]
    
    .
  




{circumflex over (q)}mi denotes the 3D unit vector {overscore (q)}/|{overscore (q)}|=(q;1)T/|(q;1)|. Similarly, Â=A/|A|. Let smi≡qmi=qm0 denote the image displacement for the m-th tracked point and Qm=(Qx, Qy, Qz)mT denote 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)T be the image coordinates of the n-the pixel position. Order the pixels from 1 to Nx, the total number of pixels. Let Ii denote the i-the intensity image and let Ini=Ii(pn) denote the image intensity at the n-the pixel position in Il. Let Pn denote the 3D point imaged at pn in the reference image, with Pn≡(Xn, Xn, Xn)T in the coordinate system of I0. Let dni denote the shift in image position from Io to Ii of 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
                    ^
                  
                  n
                  i
                
                -
                
                  
                    R
                    i
                  
                  ⁢
                  
                    q
                    m
                    o
                  
                
              
              )
            
            2
          
        
        +
        
          μ
          ⁢
          
            
              ∑
              l
            
            ⁢
            
              
                (
                
                  
                    
                      A
                      ^
                    
                    l
                    i
                  
                  -
                  
                    
                      R
                      i
                    
                    ⁢
                    
                      
                        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 are 



  
    s
    m
    i
  
  =
  
    
      
        
          Q
          
            z
            ,
            m
          
          
            -
            1
          
        
        ⁡
        
          (
          
            
              
                q
                m
              
              ⁢
              
                T
                z
                i
              
            
            -
            
              
                [
                
                  T
                  i
                
                ]
              
              2
            
          
          )
        
      
      
        1
        -
        
          
            Q
            
              z
              ,
              m
            
            
              -
              1
            
          
          ⁢
          
            T
            z
            i
          
        
      
    
    +
    
      f
      ⁡
      
        (
        
          
            R
            i
          
          ,
          
            q
            m
            i
          
        
        )
      
    
  




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-2Np translational flow vectors 



  
    Φ
    x
  
  ≡
  
    
      -
      
        [
        
          
            
              
                (
                
                  Q
                  z
                  
                    -
                    1
                  
                
                }
              
            
          
          
            
              
                {
                0
                }
              
            
          
        
        ]
      
    
    ⁢
    
       
    
    ⁢
    
      Φ
      y
    
  
  ≡
  
    
      -
      
        [
        
          
            
              
                {
                0
                }
              
            
          
          
            
              
                {
                
                  Q
                  z
                  
                    -
                    1
                  
                
                }
              
            
          
        
        ]
      
    
    ⁢
    
       
    
    ⁢
    
      Φ
      z
    
  
  ≡
  
    -
    
      [
      
        
          
            
              {
              
                
                  q
                  x
                
                ⁢
                
                  Q
                  z
                  
                    -
                    1
                  
                
              
              }
            
          
        
        
          
            
              {
              
                
                  q
                  y
                
                ⁢
                
                  Q
                  z
                  
                    -
                    1
                  
                
              
              }
            
          
        
      
      ]
    
  




and the length −2Np rotational flows 



  
    Ψ
    x
  
  =
  
    
      
        [
        
          
            
              
                {
                
                  
                    r
                    x
                    
                      (
                      1
                      )
                    
                  
                  ⁡
                  
                    (
                    q
                    )
                  
                
                }
              
            
          
          
            
              
                {
                
                  
                    r
                    y
                    
                      (
                      1
                      )
                    
                  
                  ⁡
                  
                    (
                    q
                    )
                  
                
                }
              
            
          
        
        ]
      
      ⁢
      
         
      
      ⁢
      
        Ψ
        y
      
    
    ≡
    
      
        [
        
          
            
              
                {
                
                  
                    r
                    x
                    
                      (
                      2
                      )
                    
                  
                  ⁡
                  
                    (
                    q
                    )
                  
                
                }
              
            
          
          
            
              
                {
                
                  
                    r
                    y
                    
                      (
                      2
                      )
                    
                  
                  ⁡
                  
                    (
                    q
                    )
                  
                
                }
              
            
          
        
        ]
      
      ⁢
      
         
      
      ⁢
      
        Ψ
        z
      
    
    ≡
    
      [
      
        
          
            
              {
              
                
                  r
                  x
                  
                    (
                    3
                    )
                  
                
                ⁡
                
                  (
                  q
                  )
                
              
              }
            
          
        
        
          
            
              {
              
                
                  r
                  y
                  
                    (
                    3
                    )
                  
                
                ⁡
                
                  (
                  q
                  )
                
              
              }
            
          
        
      
      ]
    
  




Then collect the shifts, smi smi into a Ni×2Np matrix 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 ΘFOV denote 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 Az more accurately than Ax, Ay by 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 using 



  
    
      A
      *
    
    =
    
      [
      
        
          
            
              R
              ⁢
              
                 
              
              ⁢
              A
            
          
        
        
          
            
              T
              ·
              A
            
          
        
      
      ]
    
  
  ,
  
    
      B
      *
    
    =
    
      
        [
        
          
            
              
                R
                ⁢
                
                   
                
                ⁢
                B
              
            
          
          
            
              
                
                  T
                  ·
                  B
                
                +
                
                  B
                  4
                
              
            
          
        
        ]
      
      .
    
  




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 by 



  
    
      
        
          A
          ′
        
        ≡
        
          
            R
            ⁡
            
              (
              
                A
                -
                
                  B
                  ⁢
                  
                    
                      T
                      ·
                      A
                    
                    
                      
                        T
                        ·
                        B
                      
                      +
                      
                        B
                        4
                      
                    
                  
                
              
              )
            
          
          .
        
      
    
    
      
        (
        3
        )
      
    
  




Now consider the line flow δAli=Ali−Al0. Define δAR=Al−R−1Ai and 



  
    
      δ
      ⁢
      
         
      
      ⁢
      
        A
        T
      
    
    ≡
    
      
        
          R
          
            -
            1
          
        
        ⁢
        
          A
          i
        
      
      -
      A
    
  
  =
  
    
      
        
          (
          
            T
            ·
            A
          
          )
        
        
          1
          -
          
            T
            ·
            B
          
        
      
      ⁢
      B
      ⁢
      
         
      
      ⁢
      so
      ⁢
      
         
      
      ⁢
      δ
      ⁢
      
         
      
      ⁢
      A
    
    =
    
      
        δ
        ⁢
        
           
        
        ⁢
        
          A
          T
        
      
      +
      
        δ
        ⁢
        
           
        
        ⁢
        
          
            A
            R
          
          .
        
      
    
  




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)·Alo exactly. We normalize all A0 to 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 δAli incorporates just two degrees of freedom. We represent δAli by 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 PUl and PLl project onto these two directions. For typical ΘFOV<1,|Al·{circumflex over (z)}|<<1 and the upper component of Al roughly 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
    
  
  ≡
  
    
      [
      
        
          
            
              {
              
                
                  A
                  x
                
                ⁢
                
                  B
                  U
                
              
              }
            
          
        
        
          
            
              {
              
                
                  A
                  x
                
                ⁢
                
                  B
                  L
                
              
              }
            
          
        
      
      ]
    
    ⁢
    
       
    
    ⁢
    
      Φ
      
        L
        ⁢
        
           
        
        ⁢
        y
      
    
  
  ≡
  
    
      [
      
        
          
            
              {
              
                
                  A
                  y
                
                ⁢
                
                  B
                  U
                
              
              }
            
          
        
        
          
            
              {
              
                
                  A
                  y
                
                ⁢
                
                  B
                  L
                
              
              }
            
          
        
      
      ]
    
    ⁢
    
       
    
    ⁢
    
      Φ
      z
    
  
  ≡
  
    [
    
      
        
          
            {
            
              
                A
                z
              
              ⁢
              
                B
                U
              
            
            }
          
        
      
      
        
          
            {
            
              
                A
                z
              
              ⁢
              
                B
                L
              
            
            }
          
        
      
    
    ]
  




where these are length-2NL vectors. BU≡B·PU and BL≡B·PL are the upper and lower components of B. Similarly, define the rotational flows 



  
    
      
        
          Ψ
          
            L
            ⁢
            
               
            
            ⁢
            x
          
        
        ≡
        
          
            [
            
              
                
                  
                    {
                    
                      
                        P
                        U
                      
                      ·
                      
                        (
                        
                          x
                          ×
                          A
                        
                        )
                      
                    
                    }
                  
                
              
              
                
                  
                    {
                    
                      
                        P
                        L
                      
                      ·
                      
                        (
                        
                          x
                          ×
                          A
                        
                        )
                      
                    
                    }
                  
                
              
            
            ]
          
          ⁢
          
             
          
          ⁢
          
            Ψ
            
              L
              ⁢
              
                 
              
              ⁢
              y
            
          
        
        ≡
        
          
            [
            
              
                
                  
                    {
                    
                      
                        P
                        U
                      
                      ·
                      
                        (
                        
                          y
                          ×
                          A
                        
                        )
                      
                    
                    }
                  
                
              
              
                
                  
                    {
                    
                      
                        P
                        L
                      
                      ·
                      
                        (
                        
                          y
                          ×
                          A
                        
                        )
                      
                    
                    }
                  
                
              
            
            ]
          
          ⁢
          
             
          
          ⁢
          
            Ψ
            z
          
        
        ≡
        
          [
          
            
              
                
                  {
                  
                    
                      P
                      U
                    
                    ·
                    
                      (
                      
                        z
                        ×
                        A
                      
                      )
                    
                  
                  }
                
              
            
            
              
                
                  {
                  
                    
                      P
                      L
                    
                    ·
                    
                      (
                      
                        z
                        ×
                        A
                      
                      )
                    
                  
                  }
                
              
            
          
          ]
        
      
    
    
      
        (
        5
        )
      
    
  




Let Λ be the N1×2NL matrix 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 smi for 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)nT represent the gradient of the image intensity for I0, with some appropriate definition for a discrete grid of pixels. Let ΔIni define 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) 

This and the previous equation imply that −ΔIni=∇In·dni≈Zn−1(∇In·pnTzi−∇In·[Ti]2)+∇In·(ωxirn(1)+ωyir(2)+ωyirn(3)). Define the three length −Nx translational flow vectors Φlx≡−{Z−1Ix}, Φly≡−{Z−1Iy}, Φlz≡{Z−1(∇I·p)}, and the length −Nx rotational flows, ωlx≡−{∇I·r(1)(p)}, ωly≡−{∇I·r(2)(p)}, ωlz≡−{∇I·r(3)(p)}. Let Δ be a N1×Nx matrix, where each row corresponds to an image i and equals {ΔIi}T. Then

Δ≈{Tx}ΦlxT+{Ty}ΦlyT+{Tz}ΦlzT+{ωz}ΨlxT+{ωy}ΨlyT+{Ψz}ΨlzT.  (8) 

Reconstruction

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 ωL and ω1 are 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)×Ntot matrix H annihilating the three {overscore (ω)}x,y,z. Computing and multiplying by H cost O(Ntot) computation. Define the N1×Ntot matrix {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. Define 



  
    D
    _
  
  ≡
  
    
      C
      
        -
        
          1
          2
        
      
    
    ⁢
    
      D
      _
    
    ⁢
    
       
    
    ⁢
    
      
        H
        T
      
      .
    
  




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 that 



  
    
      
        
          
            D
            _
          
          CH
        
        ≈
        
          
            
              C
              
                -
                
                  1
                  2
                
              
            
            ⁢
            
              {
              
                T
                x
              
              }
            
            ⁢
            
              
                Φ
                _
              
              T
            
            ⁢
            
              H
              T
            
          
          +
          
            
              C
              
                -
                
                  1
                  2
                
              
            
            ⁢
            
              {
              
                T
                y
              
              }
            
            ⁢
            
              
                Φ
                _
              
              y
              T
            
            ⁢
            
              H
              T
            
          
          +
          
            
              C
              
                -
                
                  1
                  2
                
              
            
            ⁢
            
              {
              
                T
                z
              
              }
            
            ⁢
            
              
                Φ
                _
              
              z
              T
            
            ⁢
            
              
                H
                T
              
              .
            
          
        
      
    
    
      
        (
        9
        )
      
    
  




    {overscore (D)}CH is 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)T where 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, Bz and Z−1 from the {overscore (Φ)}α in this system of equations to get 3Ntot linear 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, Bz and Z−1 from [{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] and 



  
    
      D
      _
    
    CH
  
  ≈
  
    
      
        C
        
          -
          
            1
            2
          
        
      
      ⁢
      
        {
        
          T
          x
        
        }
      
      ⁢
      
        
          Φ
          _
        
        T
      
      ⁢
      
        H
        T
      
    
    +
    
      
        C
        
          -
          
            1
            2
          
        
      
      ⁢
      
        {
        
          T
          y
        
        }
      
      ⁢
      
        
          Φ
          _
        
        y
        T
      
      ⁢
      
        H
        T
      
    
    +
    
      
        C
        
          -
          
            1
            2
          
        
      
      ⁢
      
        {
        
          T
          z
        
        }
      
      ⁢
      
        
          Φ
          _
        
        z
        T
      
      ⁢
      
        H
        T
      
    
  




to recover the translations.

        5) We recover the rotations ωxi, ωyi, ωzi from 



  
    
      
        ω
        x
        i
      
      ⁢
      
        
          Ψ
          _
        
        xn
      
    
    +
    
      
        ω
        yx
        i
      
      ⁢
      
        
          Ψ
          _
        
        yn
      
    
    +
    
      
        ω
        zx
      
      ⁢
      
        
          Ψ
          _
        
        zn
      
    
  
  =
  
    
      
        C
        
          -
          
            1
            2
          
        
      
      ⁢
      
        
          D
          _
        
        n
        i
      
    
    -
    
      
        (
        
          
            C
            
              -
              
                1
                2
              
            
          
          ⁡
          
            (
            
              
                
                  {
                  
                    T
                    x
                  
                  }
                
                ⁢
                
                  
                    Φ
                    _
                  
                  T
                
              
              +
              
                
                  {
                  
                    T
                    y
                  
                  }
                
                ⁢
                
                  
                    Φ
                    _
                  
                  y
                  T
                
              
              +
              
                
                  {
                  
                    T
                    z
                  
                  }
                
                ⁢
                
                  
                    Φ
                    _
                  
                  z
                  T
                
              
            
            )
          
        
        )
      
      n
      i
    
  




    
    


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 Steps 2-4 we weight Tzi by a factor proportional to the FOV and weight {overscore (Φ)}z by 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−1Tz and 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 dni by 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 Tz are 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
    
    ⁢
    
      
        
          Φ
          _
        
        La
        T
      
      ⁢
      
        H
        T
      
      ⁢
      
        P
        s
      
      ⁢
      H
      ⁢
      
         
      
      ⁢
      
        
          Φ
          _
        
        la
      
    
  
  ,




will also produce small overlaps 



  
    ∑
    a
  
  ⁢
  
    
      
        Φ
        _
      
      La
      T
    
    ⁢
    
      H
      T
    
    ⁢
    
      P
      s
    
    ⁢
    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 that 



  
    
      
        B
        T
      
      ⁢
      
        H
        B
      
      ⁢
      B
    
    =
    
      
        ∑
        a
      
      ⁢
      
        
          
            Φ
            _
          
          La
          T
        
        ⁢
        
          H
          T
        
        ⁢
        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−1 toward the reference image I0. In this case we then let the image displacements dni refer to the unrotated images.
        1. We then compute H and 



  
    
      Δ
      CH
    
    ≡
    
      
        C
        
          -
          
            1
            2
          
        
      
      ⁢
      Δ
      ⁢
      
         
      
      ⁢
      
        H
        T
      
    
  
  ,




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−1 from the {overscore (Φ)}α in this system of equations to get 3Np linear 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−1 from (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
        
          -
          
            1
            2
          
        
      
    
    ⁢
    W
    ⁢
    
       
    
    ⁢
    
      Ψ
      T
    
  
  =
  
    -
    
      
        
          C
          
            -
            
              1
              2
            
          
        
        ⁡
        
          (
          
            
              
                T
                _
              
              ⁢
              
                 
              
              ⁢
              
                Φ
                T
              
            
            +
            Δ
          
          )
        
      
      .
    
  




Algorithm Analysis

    
    


Step 0.
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
          
            -
            
              1
              2
            
          
        
      
      ⁢
      Δ
    
    ⁢
    
      ≈
      
        -
        
          1
          2
        
      
    
    ⁢
    
      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
    
      -
      
        1
        2
      
    
  




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.
Step 2a.
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.


First we define {overscore (P)}≡HTH, a projection matrix. Multiplying (7) by HT gives HTS(3)U≈{overscore (P)}Φ, which implies

Φ≈HTS(3)+ΨΩ

We eliminate the Zn−1 from (8) above and solve directly for U and Ω. Denote the columns of U by U≡[U1 U2 U3], and similarly for Ω. Let U3′≡s−1U3 and Ω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′]n and 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−1 from 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−1 multi-frame estimates. To get a better multi-frame estimate, one can first re-compute I0 and ∇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
      
        -
        
          1
          2
        
      
    
    ⁢
    Δ
    ⁢
    
       
    
    ⁢
    
      
        H
        T
      
      .
    
  




Up to unknown rotational flows, ΔCH(3)H gives the intensity shifts due to the translational displacements: 



  
    
      C
      
        -
        
          1
          2
        
      
    
    ⁢
    Δ
  
  ≈
  
    
      
        Δ
        CH
        
          (
          3
          )
        
      
      ⁢
      H
    
    +
    
      W
      ⁢
      
         
      
      ⁢
      
        Ψ
        T
      
    
  




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



  
    
      
        C
        
          -
          
            1
            2
          
        
      
      ⁢
      
        T
        _
      
    
    ≈
    
      
        M
        
          (
          3
          )
        
      
      ⁢
      
        U
        T
      
    
  
  ,




where UT is also a 3×3 matrix. The algorithm recovers UT and {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−1 one must replace Ψ in (8) by [Ψ, ΨF], and Ω becomes a 4×3 matrix. The modified Step 2a can recover the Z−1 except 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-Np vectors, 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 Mi is a 3×3 matrix (a homography), Sm is the structure 3-vector giving the position of the m-th point on the plane (in some arbitrary coordinate system), and λmi is 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)}mi into a single 3N1×Np matrix Γ. 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 λmi minimizing the |Γ−M(3)S(3)|. Return to Step 2.
    
    


After convergence, M(3) contains the Mi estimates. 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 Mi if 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 Ki be 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′)−T RK≡MH and T→KT, where K and K′ are the calibration matrices for the first and second image. MH is a general homography and is the inverse transpose of the analogous homography for points. The first order approximation of the new version of 



  
    A
    ′
  
  =
  
    R
    ⁡
    
      (
      
        A
        -
        
          B
          ⁢
          
            
              T
              ·
              A
            
            
              
                T
                ·
                B
              
              +
              
                B
                4
              
            
          
        
      
      )
    
  




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−T is 1−δMHT, one can easily define H to annihilate the first-order flows for both points and lines. Represent 



  
    
      M
      H
    
    ≡
    
      [
      
        
          
            F
          
          
            G
          
        
        
          
            
              J
              T
            
          
          
            0
          
        
      
      ]
    
  
  ,




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 δMH are given in the following table:





















F1,1
F1,2
F2,2
F2,2
G1
G2
J1
J2
























Ax:
Ax
Ay
0
0
Az
0
0
0


Ay:
0
0
Ax
Ay
0
Az
0
0


Az:
0
0
0
0
0
0
Ax
Ay


qx:
−x
0
−y
0
−x2
−xy
−1
0


qy:
0
−x
0
−y
−xy
−y2
0
−1









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−1 where 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−1 given 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 to FIG. 1, which illustrates in block-diagram form a computer hardware system incorporating the invention. As indicated therein, the system includes a video source 101, whose output is digitized into a pixel map by a digitizer 102. The digitized video frames are then sent in electronic form via a system bus 103 to a storage device 104 for access by the main system memory during usage. During usage the operation of the system is controlled by a central-processing unit, (CPU) 105 which 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, 106 and graphics board 107.
The user interacts with the system by way of a keyboard 108 and or a mouse 109 or 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 to FIG. 2 the 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 step 201. Me images are then digitized 202 for 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 determined 203; 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 constructed 204 and 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 calculated 205. Recovering the 3D structure from the 3D structure matrix using SVD by solving a linear equation 206. The camera motion can then be recovered using the recovered 3D structure 207.
While there has been shown and described what is considered to be preferred embodiments of the invention, it will, of course, be understood that various modifications and changes in the form or detail could readily be made without departing from the spirit of the invention. It is therefore intended that the invention be not limited to the exact forms described and illustrated, but should be constructed to cover all modifications that may fall within the scope of the appended claims.

## Claims

1. An image processing method for recovery of a scene structure from successive image data, the method comprising the steps of:
(a) computing rotational motion in the successive image data using a combination of first rotational flow vectors derived from a set of tracked points collected from the successive image data second rotational flow vectors derived from a sat of tracked lines collected from the successive image data, where the first rotational flow vectors are represented by 



  
    
      Ψ
      x
    
    ≡
    
      [
      
        
          
            
              {
              
                
                  r
                  x
                  
                    (
                    1
                    )
                  
                
                ⁡
                
                  (
                  q
                  )
                
              
              }
            
          
        
        
          
            
              {
              
                
                  r
                  y
                  
                    (
                    1
                    )
                  
                
                ⁡
                
                  (
                  q
                  )
                
              
              }
            
          
        
      
      ]
    
  
  ,
  
    
      Ψ
      y
    
    ≡
    
      [
      
        
          
            
              {
              
                
                  r
                  x
                  
                    (
                    2
                    )
                  
                
                ⁡
                
                  (
                  q
                  )
                
              
              }
            
          
        
        
          
            
              {
              
                
                  r
                  y
                  
                    (
                    2
                    )
                  
                
                ⁡
                
                  (
                  q
                  )
                
              
              }
            
          
        
      
      ]
    
  
  ,
  
    
      Ψ
      z
    
    ≡
    
      [
      
        
          
            
              {
              
                
                  r
                  x
                  
                    (
                    3
                    )
                  
                
                ⁡
                
                  (
                  q
                  )
                
              
              }
            
          
        
        
          
            
              {
              
                
                  r
                  y
                  
                    (
                    3
                    )
                  
                
                ⁡
                
                  (
                  q
                  )
                
              
              }
            
          
        
      
      ]
    
  





where r(1), r(2), r(3) are three point rotational flows with respect to an image position q of one of the tracked points and where the second rotational flow vectors derived from the set of tracked lines are represented by projections in two directions that take into account differences in noise of a measured line in different directions and where the second rotational flow vectors are represented by 



  
    
      Ψ
      Lx
    
    ≡
    
      [
      
        
          
            
              {
              
                
                  P
                  U
                
                ·
                
                  (
                  
                    
                      x
                      ^
                    
                    ×
                    A
                  
                  )
                
              
              }
            
          
        
        
          
            
              {
              
                
                  P
                  L
                
                ·
                
                  (
                  
                    
                      x
                      ^
                    
                    ×
                    A
                  
                  )
                
              
              }
            
          
        
      
      ]
    
  
  ,
  
    
      Ψ
      Ly
    
    ≡
    
      [
      
        
          
            
              {
              
                
                  P
                  U
                
                ·
                
                  (
                  
                    
                      y
                      ^
                    
                    ×
                    A
                  
                  )
                
              
              }
            
          
        
        
          
            
              {
              
                
                  P
                  L
                
                ·
                
                  (
                  
                    
                      y
                      ^
                    
                    ×
                    A
                  
                  )
                
              
              }
            
          
        
      
      ]
    
  
  ,
  
    
      Ψ
      Lz
    
    ≡
    
      [
      
        
          
            
              {
              
                
                  P
                  U
                
                ·
                
                  (
                  
                    
                      z
                      ^
                    
                    ×
                    A
                  
                  )
                
              
              }
            
          
        
        
          
            
              {
              
                
                  P
                  L
                
                ·
                
                  (
                  
                    
                      z
                      ^
                    
                    ×
                    A
                  
                  )
                
              
              }
            
          
        
      
      ]
    
  




where A represents one of the tracked lines and PU and PL represent unit vectors projecting on the two directions;
(b) constructing a shift data matrix for the tracked lines and the tracked points that compensates for the rotational motion in the successive image data; and 
(c) decomposing the shi data matrix into a translational motion matrix and a structure matrix and recovering the scene structure by solving for the structure matrix.

2. The image processing method of claim 1 wherein ith rotational motion is compensated through linear elimination of rotational dependence in the shift data matrix.

3. The image processing method of claim 1 wherein the shift data matrix is decomposed using singular value docomposition.

4. The image processing method of claim 1 wherein components of the shift data matrix ae weighted to account for greater accuracy of measurement of the components.

5. The image processing method of claim 1 wherein the method is iterated until it converges to a reconstruction of the scene structure.

6. An image processing method for recovery of a scene structure from successive image data, the method comprising the steps of:
(a) parameterizing a set of tracked lines collected from the successive image data; 
(b) computing rotational motion in the successive image data using rotational flow vectors derived from the tracked lines where the rotational flow vectors am represented by projections in two directions that take into account differences in noise of a measured line in different directions and where the rotational flow vectors are represented by 



  
    
      Ψ
      Lx
    
    ≡
    
      [
      
        
          
            
              {
              
                
                  P
                  U
                
                ·
                
                  (
                  
                    
                      x
                      ^
                    
                    ×
                    A
                  
                  )
                
              
              }
            
          
        
        
          
            
              {
              
                
                  P
                  L
                
                ·
                
                  (
                  
                    
                      x
                      ^
                    
                    ×
                    A
                  
                  )
                
              
              }
            
          
        
      
      ]
    
  
  ,
  
    
      Ψ
      Ly
    
    ≡
    
      [
      
        
          
            
              {
              
                
                  P
                  U
                
                ·
                
                  (
                  
                    
                      y
                      ^
                    
                    ×
                    A
                  
                  )
                
              
              }
            
          
        
        
          
            
              {
              
                
                  P
                  L
                
                ·
                
                  (
                  
                    
                      y
                      ^
                    
                    ×
                    A
                  
                  )
                
              
              }
            
          
        
      
      ]
    
  
  ,
  
    
      Ψ
      Lz
    
    ≡
    
      [
      
        
          
            
              {
              
                
                  P
                  U
                
                ·
                
                  (
                  
                    
                      z
                      ^
                    
                    ×
                    A
                  
                  )
                
              
              }
            
          
        
        
          
            
              {
              
                
                  P
                  L
                
                ·
                
                  (
                  
                    
                      z
                      ^
                    
                    ×
                    A
                  
                  )
                
              
              }
            
          
        
      
      ]
    
  





where A represents one of the tracked lines and PU and PL represent unit vectors projecting on the two directions;
(c) constructing a shift data matrix for the tracked lines that compensates for the rotational motion in the successive image data; and 
(d) decomposing the shift data matrix into a translational motion matrix and a structure matrix and recovering the scene structure by solving for the structure matrix.

7. The image processing method of claim 6 wherein the rotational motion is compensated through linear elimination of rotational dependence in the shift data matrix.

8. The image method of claim 6 wherein the shift data matrix is decomposed using singular value decomposition.

9. The image method of claim 6 wherein components of the shift data matrix are weighted to account for greater accuracy of measurement of the components.

10. The image processing method of claim 6 wherein the method is iterated until it converges to a reconstruction of the scene structure.