Patent Publication Number: US-9898822-B2

Title: Method and system for reconstructing scenes as 3D models from sequences of images using constraint satisfaction

Description:
FIELD OF THE INVENTION 
     The invention is generally related to computer vision, and more particularly to reconstructing scenes from sequences of images using constraint satisfaction to determine relative poses and to reconstruct 3D points and 3D lines. 
     BACKGROUND OF THE INVENTION 
     In computer vision applications, estimation of geometrical characteristics of scenes and objects, such as relative pose, absolute pose and three-dimensional (3D) models, rely on underlying constraints. For example, graph-rigidity based on inter-point distances can be used for 3D reconstruction. Cheirality is one of the key concepts that enables one to identify a correct reconstruction from ambiguous ones arising from mirror reflection. In simple words, the cheirality constraint enforces that the scene is in front of the camera rather than behind it. 
     Divide and Concur (D&amp;C) is a constraint satisfaction method based on difference-map (DM) dynamics. The D&amp;C method is closely related to iterative optimization methods. In particular, the D&amp;C method can be characterized as a message-passing method and it is similar to Belief Propagation (BP). The D&amp;C method can be characterized as a special case of an alternating direction method of multipliers (ADMM) method. Furthermore, the D&amp;C method is a generalization of iterative projection methods such as Douglas-Rachford for convex optimization, and Fienup&#39;s Hybrid input output (HIO) method for phase retrieval. 
     Generally, 3D modeling methods use the following steps. First, one estimates the poses of the cameras, and then outliers are removed. To this end, minimal approaches are critical. After obtaining some initial solution for camera poses and 3D points, bundle-adjustment can be applied. 
     In a 5-point solver, several challenging minimal problems are addressed using algebraic solvers. Algebraic solvers are not guaranteed to solve any polynomial system. For example, certain minimal problems such as 3-view 4-point are notoriously difficult to solve. There are other challenging problems, such as estimating the poses from six line correspondences in three views, and minimal methods for general cameras. 
     General cameras refer to a broad class of sensors where the projection rays do not necessarily pass through an optical center of the camera. This generality creates greater challenges when one tries to solve minimal problems. For example, the relative pose between two general cameras produces 64 solutions after solving a large polynomial system. 
     SUMMARY OF THE INVENTION 
     The embodiments of invention provide a method and system for reconstructing a scene as a three-dimensional (3D) model from a sequence of images using constraint satisfaction to determine relative poses and to reconstruct 3D points and 3D lines. The constraint satisfaction is based on a difference-map (DM) search strategy. The method does not require any initialization and can determine both discrete and continuous variables. Existing approaches require initialization and are usually customized for either discrete or continuous variables, but not both. 
     The method can be considered as a special case of an alternating direction method of multipliers (ADMM). The embodiments model computer vision problems, such as a line-labeling problem, 3D reconstruction using point or line correspondences for both perspective, and general cameras, line-based 3D reconstruction using geometrical constraints. 
     In contrast to other methods that use a single image or fewer constraints, the embodiments can concurrently handle all the available constraints and images. 
     Iterative methods are not new in computer vision. In fact, one of the most applied methods in geometric computer vision is bundle-adjustment. Despite its significant impact on handling city-scale 3D models, local minima are still a serious problem for bundle adjustment. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  is a block diagram of a method and system for reconstructing a scene as a three-dimensional (3D) model from a sequence of images according to embodiments of the invention; 
         FIGS. 2A, 2B and 2C  are line drawings with labels according to embodiments of the invention; 
         FIGS. 3A, 3B, 3C and 3D  are schematics of T, W, Y and L labeled junctions according to embodiments of the invention; 
         FIG. 4  schematically shows 3D reconstruction from multiple perspective images; 
         FIG. 5  schematically shows 3D reconstruction from multiple perspective images for general cameras; 
         FIG. 6  schematically shows 3D reconstruction with Cheirality and rigidity constraints; and 
         FIG. 7 . schematically shows 3D reconstruction of lines from multiple images using plane constraints on points lying on intersecting lines. 
     
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
     As shown in  FIG. 1 , the embodiments of our invention provide a method and system for reconstructing a scene as a three-dimensional (3D) model from a sequence of images using constraint satisfaction to determine relative poses and the reconstruction of 3D points and lines. 
     A sequence of images  102  is acquired of a scene  101  with a camera. Feature point matches or line matches are extracted  110  from the sequence of images using, for example, a conventional feature matching procedure, such as a scale-invariant feature transform (SIFT). 
     Next, variables are initialized  120 , e.g., camera optical centers and 3D points, using random values. We make  130  n copies of the variables depending on the availability of n constraints  130 . We project  140  the copies to satisfy each of the constraints. The copies of the same variables are forced to agree using a concur constraint by averaging the replicas. We iterate the steps  140  and  150  until convergence. The average values of the variables after convergence provide the 3D model  160 . 
     The embodiments of our invention apply to 3D modeling that uses both discrete and continuous optimization techniques, and Divide and Concur (D&amp;C) for constraint satisfaction. In particular, belief propagation is used for solving rotation matrices, and Levenberg Marquardt iteration is used for continuous optimization. 
     The embodiments provide solutions for both discrete (line-labeling problem), and continuous optimization (3D modeling). Conventional iterative methods, such as the Douglas-Rachford method, are generally applied on convex problems. The success of difference-map on phase-retrieval problem indicates that these methods can also be useful for non-convex problems. 
     The Difference-Map Dynamics 
     In the D&amp;C method, each variable is replicated as many as the number of constraints that uses it. Rather than working directly with the original variables of the constraint satisfaction problem, the DC method tries to determine a setting of the replica values that satisfies all the original constraints, with additional constraints that all values of the replicas of the same variables must be equal. The original ones are called Divide constraints, while the equalities are called Concur constraints. The DC method operates by alternating the Divide and Concur constraints. 
     The D&amp;C projections are combined according to Difference-Map dynamics. We briefly describe how the difference-map works and how it performs better than a simple alternating projection method. The difference-map overshoots the replicas after satisfying the Divide constraints. Then, the overshoots are subtracted from the replicas. This allows the D&amp;C method to avoid the traps (fixed points that are not solutions to the problem), which causes serious issues for the simple alternating projection method. 
     Each constraint on a point x, i.e., a vector of variables, can be expressed as a set S and the constraint satisfaction problem is equivalent to finding a point in the set x′εS such that p′=argmin p , |p−p′|, where |p−p′| is the Euclidean distance. The associated operator P: x→x′, {x′εS} is referred to as the projection. In other words, the projection operator moves the variables by a minimal distance so that the variables satisfy a particular constraint. 
     Let us consider a set of N constraints, each denoted by a set S i , i={1, 2, . . . , n}. Our goal is to determine a point {circumflex over (p)} that lies in the intersection of all these sets. In other words, {circumflex over (p)}εS i , ∀i={1, 2, . . . , n}. We first make n replicas of the variables at a time t as
 
 y   t =( x   t:1   ,x   t:2   , . . . ,x   t:n ).  (1)
 
     Let us denote the projections for each of these constraints as P i , i={1, . . . , n}. We use P D  and P C  to denote “Divide” and “Concur” projections respectively. The divide projection P D  moves the individual replicas to satisfy each of the individual constraints. In other words, the Divide projection refers to the simultaneous application of each of the individual constraints P i  on each of the individual replicas x t:i  as
 
 P   D ( y   t )=( P   1 ( x   t:1 ), P   2 ( x   t:2 ), . . . , P   N ( x   t:n )).  (2)
 
     After the application of the constraints, we make the replicas agree with each other. This is done by the “Concur” projection
 
 P   C ( P   D ( y   t ))=(   t , t , . . . , t   ).  (3)
 
where
 
               x   _     =           ∑     i   =   1     n     ⁢       P   i     ⁡     (     x     t   :   i       )         n     .           
In the next iteration, we have the following replicas:
 
 y   t+1 =(   x , x , . . . , x   ).  (4)
 
     The above method can be written in a single equation:
 
 y   t+1   =P   C ( P   D ( y   t )).  (5)
 
     This is generally referred to as alternating projections method. This method can be used in convex problems, whereas in non-convex problems, the method is prone to getting struck in fixed points which do not correspond to the solutions of the problem. We refer to a point y* as a fixed point when y t+1 =y t =y*. The alternating projections method using the DM can handle non-convex constraints in diffraction imaging problems
 
 y   t+1   =P   C ( y   t   over )−( P   D ( y   t )− y   t ).  (6)
 
     The above difference-map dynamics is a three-step process:
         1. We move the replicas by twice the amount given by the Divide projection. We refer to these as overshoot values y t   over =y t +2(P D (y t )−y t ).   2. We apply the Concur projection on the overshoot values by y t   conc =P c (y t   over ).   3. Then, the overshoot is removed from the concurred values to obtain the replicas for the next iteration y t+1 =y t   conc −(P D (y t )−y t ).       

     In contrast to the alternating projection method, if a fixed point y* given by y t+1 =y t =y* is reached by the DM dynamics, then this fixed point corresponds to a solution for the problem. It is important to note that the fixed points are not the solution to the problem. In general, it is not a single point. Instead, there is a continuum of fixed points associated with every solution. In order to determine the solution from the fixed points, we perform the Divide projection one more time as y sol =P D (y*). 
     Modeling Vision Problems as CSPs: 
     The D&amp;C method can be used to solve several problems in geometric computer vision. In order to apply D&amp;C, we first identify the constraints in the problem, and a method to enforce each constraint individually. 
     Conventional Line-Labeling Problem 
     Given a line drawing composed of line segments and junctions, our goal is to identify the junctions and label the line segments. 
       FIGS. 2A, 2B and 2C  show schematic of line drawings with four labels (+, −, →, ←), where + and − represent convex (+), concave (−) line segments, and an occluding line segment with an arrow directed such that the occluding face is on the right and the occluded face is on the left (←,→). Each line segment in a line drawing corresponds to an edge between two faces in the three-dimensional object that the line drawing represents, or is a boundary between a face and the background. 
     Junction Types 
     Every junction is an intersection of two or more edges. The junction type can be inferred easily using the number and orientations of the edges incident to the junction. We assume that every junction type is identified prior to the labeling. We consider the following junctions L, Y, T, W, M, K, X for trihedral and tetrahedral objects. An object is trihedral if three faces meet at every junction. The object is an extended trihedral when exactly three planes meet at a point. 
     Junction Labeling 
     As shown in  FIGS. 3A, 3B, 3C and 3D , the labeling of a junction denotes the label assignment for the edges in a specific order, which depends on the type of junction. The junctions T, W, Y and L are shown along with the ordering of edges incident to each one of the junctions. Here, f(i)=mod(i, 3). 
     There are only a few allowed labelings for every junction type depending on the type of object to which the junction belongs. We use four alphabets (a, b, c, d) to represent the labels. These alphabets are assigned to the edges based on the incident junction. If there is an occluding edge entering a junction, then it is assigned a label a. If there is an occluding edge leaving a junction, then it is assigned a label b. Thus, we have two labels for every edge depending on the two incident junctions. The alphabets c and d refer to convex(+) and concave(−) labels. 
     A junction catalog determines what possible labels the edges of a particular junction take. For example, a Y junction in a trihedral line drawing can only take the configurations ccc, ddd, abd, bda, and dab in the order defined in  FIG. 2 . 
     For every edge, we denote the label using a state vector: (p a  p b  p c  p d ), where p a , p b , p c  and p d  are probabilities for the four possible labels and their sum is equal to 1. Let us consider a junction J with m incident edges. The junction state x t  can be represented by a 4×m vector as
 
 x   t =( p 1 a    p 1 b    p 1 c    p 1 d    p 2 a   ·pm   d ) T .  (7)
 
     Let us assume that a specific junction J has n different labelings in the junction catalog. Then, its junction catalog is given by the following matrix: 
     
       
         
           
             
               
                 
                   
                     
                       
                         x 
                         t 
                       
                       ⁡ 
                       
                         ( 
                         
                           
                             
                               1 
                             
                             
                               0 
                             
                             
                               . 
                             
                             
                               . 
                             
                             
                               . 
                             
                             
                               0 
                             
                           
                           
                             
                               0 
                             
                             
                               0 
                             
                             
                               . 
                             
                             
                               . 
                             
                             
                               . 
                             
                             
                               0 
                             
                           
                           
                             
                               0 
                             
                             
                               1 
                             
                             
                               . 
                             
                             
                               . 
                             
                             
                               . 
                             
                             
                               0 
                             
                           
                           
                             
                               1 
                             
                             
                               0 
                             
                             
                               . 
                             
                             
                               . 
                             
                             
                               . 
                             
                             
                               1 
                             
                           
                           
                             
                               . 
                             
                             
                               . 
                             
                             
                               . 
                             
                             
                               . 
                             
                             
                               . 
                             
                             
                               . 
                             
                           
                         
                         ) 
                       
                     
                     
                       4 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       m 
                       × 
                       n 
                     
                   
                   . 
                 
               
               
                 
                   ( 
                   8 
                   ) 
                 
               
             
           
         
       
     
     Note that the above matrix has only ones and zeros for the state probabilities. 
     Let L={l 1 , l 2 , . . . , l n } denote the columns of the above matrix. In the Divide projection, we have P D (x t )=argmin l     i     εY |x t −l i |. For every edge, we have two different labelings from the overshoot operations on the two junctions. 
     Let y 1   over  and y 2   over  be the two replicas after overshooting as given by: y 1   over =(p a  p b  p c  p d ) T , y 2   over =(q a  q b  q c  q d ) T . Note that the labeling of a junction includes the states of all the edges incident on it. Thus y 1   over  and y 2   over  are the sub-vectors of the states of the two junctions corresponding to the chosen edge. 
     In the Concur projection, we average the two overshoots as 
     
       
         
           
             
               
                 
                   
                     
                       P 
                       C 
                     
                     ⁡ 
                     
                       ( 
                       
                         y 
                         over 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       
                         ( 
                         
                           
                             
                               
                                 
                                   
                                     p 
                                     a 
                                   
                                   + 
                                   
                                     q 
                                     a 
                                   
                                 
                                 2 
                               
                             
                             
                               
                                 
                                   
                                     p 
                                     b 
                                   
                                   + 
                                   
                                     q 
                                     b 
                                   
                                 
                                 2 
                               
                             
                             
                               
                                 
                                   
                                     p 
                                     c 
                                   
                                   + 
                                   
                                     q 
                                     c 
                                   
                                 
                                 2 
                               
                             
                             
                               
                                 
                                   
                                     p 
                                     d 
                                   
                                   + 
                                   
                                     q 
                                     d 
                                   
                                 
                                 2 
                               
                             
                           
                         
                         ) 
                       
                       T 
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
     Because every edge is adjacent to two junctions, we have two replicas for every edge and the D&amp;C method iterates until all junctions take valid labelings from the catalog and every edge is assigned the same label by both the adjacent junctions. The states of the edges change between integers and real values during Divide and Concur projections respectively. 
     Multi-View Geometry Problems 
     We list a few basic constraints that are sufficient to solve many computer vision problems. We write each of the constraint in the form y=P(x,z), where x is the input variables, and y is the output obtained by moving x by a minimal distance to enforce the constraint. The term z refers to observed variables whose values can not be changed. 
     Rigidity y=P r (x,z): Given a set of 3D points x and a template set z, we move x to y so that y and z are related by a Euclidean transform. 
     Line constraint y=P l (x): Given a set of 3D points x, we determine y that lies on a line. 
     Plane constraint y=P π (x): Given a set of 3D points x, we determine y that lies on a plane. 
     Cheirality constraint y=P χ (x): Given three points (origin, direction vector, and a particular point), we ensure that the point is in front of the origin with respect to the direction vector. We also ensure that this point is at least one unit from the origin. 
     Directed line y=P dl (x,z): Given a set of points x and a direction vector z, we determine y that lies on a line with direction z. 
     All these constraints can be enforced. For example, in order to enforce the plane constraint y=P π (x), we first fit a plane π to the points x using least-squares and then project all the points on to this plane π to obtain y. 
     We show schematic diagrams for four geometry problems in  FIGS. 4-7 . We solve each problem by considering a set of points to denote the scene and the camera. Next, we identify a few simple constraints such as line, plane and rigidity to provide an input for the D&amp;C. We briefly describe each of these problems below. 
     n-View m-Point Perspective Case 
     In  FIG. 4 , we show the problem of 3D reconstruction from multiple perspective images. In order to solve the problem, we define the variables and the constraints. The two camera optical centers are denoted by o 1  and o 2  respectively. We consider 3D direction points d ij  at unit distances from the camera optical centers. The 3D points are denoted by P 1  and P 2  respectively. In order to solve this problem, we consider several constraints that act on these variables. We enforce a line constraint on triplets of points consisting of the camera optical center, the direction point and the associated 3D point. 
     For example, (o 1 , d 11 , p 1 ) can be one such triplet out of four. We enforce rigidity constraint for every camera. The ray directions and camera centers are rigid as per the camera calibration. During the iterations, the 3D points move around to satisfy the line constraints. We enforce rigidity constraints on the triplets (o 1 , d 11 , d 12 ) and (o2 1 , d 21 , d 22 ). The rigidity constraint allow this triplet to move around as a whole without allowing any changes in the inter-point distances among the pairs of points in the triplet. 
     On each projection ray, we enforce the Cheirality constraint that ensures that the 3D scene is reconstructed in front of the camera and also the points are at least 1 unit from the camera. The perspective reconstruction is always up to a scale. This constraint is necessary to avoid shrinking the scene to a single point. 
     n-View m-Point General Case 
     General cameras can be used to model a wide variety of camera models. In the general camera, the projection rays do not pass through a single point. Hence, we need two different 3D points to represent every ray. For example, in the first general camera in  FIG. 5 , we use o 11  and d 11  to denote a projection ray. The other constraints are similar to the ones used for perspective case. Note that there is no scale ambiguity in the case of general cameras. 
     n-View m-Line General Case 
     Given a set of line correspondences from different views of a general camera, our goal is to recover both the structure of the lines and the camera poses. In general, line-based triangulation is more challenging when compared to point-based ones. Here, we describe the problem of recovering 3D lines and camera poses for general cameras. 
     As shown in  FIG. 6 , in addition to the line through points p i,j , Cheirality and rigidity constraints, we also use additional line constraints on the 3D scene line. 
     n-View m-Line Perspective Case with Geometric Constraints 
     Given a set of line correspondences, directions in the world, and connectivity constraints, we would like to recover the structure and poses of the cameras. As shown in  FIG. 7 , we use line constraint on the points lying on the 3D world line and plane constraints on the points coining from connected line segments. 
     Although the invention has been described by way of examples of preferred embodiments, it is to be understood that various other adaptations and modifications can be made within the spirit and scope of the invention. Therefore, it is the object of the appended claims to cover all such variations and modifications as come within the true spirit and scope of the invention.