Regularizing a set of unstructured 3D points from a sequence of stereo images

Regularizing a Set of Unstructured 3D Points from a Sequence of Stereo Images Luis Alvarez, Carmelo Cuenca, and Javier S´anchez Departamento de Inform´atica y Sistemas Universidad de Las Palmas de G.C. Campus Universitario de Tafira 35017, Las Palmas {lalvarez,ccuenca,jsanchez}@dis.ulpgc.es http://serdis.dis.ulpgc.es/~{lalvarez,jsanchez} Abstract. In this paper we present a method for the regularization of a set of unstructured 3D points obtained from a sequence of stereo images. This method takes into account the information supplied by the disparity maps computed between pairs of images to constraint the regularization of the set of 3D points. We propose a model based on an energy which is composed of two terms: an attachment term that minimizes the distance from 3D points to the projective lines of camera points, and a second term that allows for the regularization of the set of 3D points by preserving discontinuities presented on the disparity maps. We embed this energy in a 2D finite element method. After minimizing, this method results in a large system of equations that can be optimized for fast computations. We derive an efficient implicit numerical scheme which reduces the number of calculations and memory allocations. 1 Introduction This paper deals with the problem of 3D geometry reconstruction from multiple 2D views. Recently, a new accurate technique based on a variational approach has been proposed in [7],[8]. Using a level set approach, this technique optimizes a 3D surface by minimizing an energy that takes into account the surface regularity as well as the projection of the surface on different images. In this paper we propose a different approach which is also based on a variational formulation but only using a disparity estimation between images and without defining explicitly any 3D surface. We will assume that the cameras are calibrated in the strong sense (see [6], [9] or [11] for more details). In the last years, very accurate techniques to estimate the disparity map in a stereo pair of images have been proposed. To extend these techniques to the case of multiple views is not a trivial problem. Roughly speaking the 3D geometry estimation that we propose can be divided in the following steps: –For each pair of consecutive images, we estimate a dense disparity map using the accurate technique developed in [1]. We estimate such disparity 2 Luis Alvarez, Carmelo Cuenca and Javier S´anchez map forward and backward, that is, from one image to the next one and in the opposite direction. –We estimate sequences of corresponding points across the multiple view image sequence. Basically, we try to connect points between images following the disparity map estimations. We select sequences of correspondent points for which the forward and backward disparity estimations are coherent. –From each selected corresponding points sequence we recover a 3D point by intersecting the projection lines of the points in the sequence. By collecting the 3D points obtained from each sequence we recover an unstructured set of 3D points. –Typically, the recovered set of 3D points is noisy, because of errors in the camera calibration process, errors in the disparity estimations, errors in the corresponding point sequences computations, etc., so some kind of regularization is needed. In this paper, we propose a new variational model to smooth the unstructured set of 3D points. This regularization model is based on the 2D image information and does not require to define any kind of geometric relation between the 3D points. The proposed technique provides a smooth set of unstructured 3D points. In this paper we do not address the problem of defining one or several surfaces fitting the set of points. Such surfaces could be recovered using some standard methods like Alpha shapes [5], Ball pivot [4] or Voronoi filtering [3]. We notice that most of such techniques require the collection of points to be smooth enough to recover the surface. So the regularization step we propose is necessary to improve the results of such techniques. The regularization model we propose, which is the main contribution of the paper is based on a variational approach. This model is designed in order to maintain the final 3D regularized surface next to the original surface and also to enable a regularization by preserving discontinuities on the disparity maps. The regularization is carried out by means of an operator which is similar to the Nagel–Enkelmann’s operator [12]. This operator has already proven its efficiency in other fields like stereoscopic reconstruction [1], optical flow estimation [2], etc. We have modified this operator to include the information given by the forward and backward optical flows of every camera. In our case it is convenient to preserve these discontinuities because disparity maps represent the depth variation of the set of 3D points. For every camera we will have a set of 2D projection points (extracted from the selected sequences of points described above) forming a grid. These grids are not necessary square and equally distributed as we would expect for regular images. On the contrary we will have in most cases some distributed points on the cameras with float precision. For this, we embed the energy in a finite element method and express the set of 3D points in terms of some basis functions. We will use a Delaunay triangulation for every camera and every point on a mesh will have a correspondent 3D point. Deriving this energy and searching for the minimum yields a system of as many equations as 3D points are there in the set. The system matrix is a sparse Lecture Notes in Computer Science 3 Fig. 1. Notation matrix that will only have no-null values in some special locations. To solve the system we first arrange the matrices to optimize their size in memory and then implement a Gauss-Seidel numerical scheme that is very efficient in the number of calculations. Finally we present two experimental results to evaluate the method. In the first experiment we use a single stereo pair of a face. We show the stereo pair and the reconstruction obtained from the method explained in [1]. Then we show the Delaunay triangulations of both cameras and finally some results of applying our regularizing method with different values for the parameters. In the other experiment we use larger sequences of stereo images. We show different results from several points of view. In Sect. 2 the notation is explained. In Sect. 3 we study the model proposing an energy with the attachment and regularizing term and explaining the operator derived from the Nagel–Enkelmann operator. In Sect. 4 we explain how to embed the method into the 2D finite element approach. In Sect. 5 we present the numerical experiments. 2 Notation As we will see later this method involves complex notation that is mainly affected by three different aspects: we use some notation for the camera configuration, notation for the stereo matching process and notation for the representation of the 3D points. We will focus on Fig. (1) to explain all these aspects. We start explaining the camera notation. We suppose there is a set of projective cameras, ci(i= 1, .., Nc), Ncdenoting the number of cameras. The projec- 4 Luis Alvarez, Carmelo Cuenca and Javier S´anchez tion model assumed for the cameras is the projective one. This model projects a 3D point X= [(X0, X1, X2)]tto a 2D image point xc= [xc, yc]tthrough a 3×4 projection matrix Pcvia sˆ xc=Pcˆ X, where sis a nonzero scale factor and the notation ˆ pis such that if p= [x1, x2, . . . , xn]tthen ˆ p= [x1, x2, . . . , xn,1]t. We note by {Xj}j=0,...,Np−1the unstructured set of 3D points, where Np is the number of points considered. Each point Xjis projected in the camera system following a sequence of corresponding points. We note by (xci(Xj), yci(Xj)) i= 0, ..., N(Xj)−1 (1) the sequence of corresponding points associated to the 3D point Xj.N(Xj) represents the number of cameras where Xjis projected. We also note by Np(c) the number of selected sequences of corresponding points passing through camera c. We may compute from the projection matrices Pc(see [6]) the focuses Fc= (Fc,0, F c,1, F c,2) for every camera c. For every 2D projection xcand for every camera cwe have an unitary vector ¯ lc(x, y)=(lc,0(x, y), lc,1(x, y), lc,2(x, y)) in the direction of focus Fc.We denote by ¯ Rc(x, y) the line passing through focus Fcand having ¯ lc(x, y) as directional vector. Let us now talk about the notation for the stereo aspects. For every pair of consecutive cameras on the sequence we have the estimated optical flow from camera cinto the previous camera c−1, ¯ hc −(x, y) = (uc −(x, y), vc −(x, y)), and the corresponding flow from camera cinto the following, ¯ hc +(x, y)=(uc +(x, y), vc +(x, y)). For every camera cwe define a function ¯ Xc:R2−→ R3that puts in correspondence bi-dimensional points (x, y) on the projection plane of camera cinto 3Dpoints (Xc,0(x, y), Xc,1(x, y), Xc,2(x, y)). Let us call X= (X0, X1, X2) the 3Dpoint obtained by the minimum distance to the set of lines ¯ Rc(x, y) that passes through focus Fcand has ¯ lc(x, y) as directional vector for the cameras in where there is a correspondent projection. Suppose that every 3Dpoint is visible from N(X) cameras. These cameras are denoted by ci(X) for i= 0, .., N(X)−1. The pixel coordinates for Xon camera ciis denoted by (xci(X), yci(X)).Therefore we have X=¯ Xc(xci(X), yci(X)) i= 0, ..., N(X)−1.(2) Finally, the derivatives of a function fare represented as df dx =fxand df dy =fy and ∇Xis the gradient of function X. 3 The 3D Regularizing Model We tackle the problem of 3D regularizing by means of a variational approach in where the solution for the minimization of a global energy is the regularized set of 3D points. Our method regularizes a set of 3D points by constraining the process with the information given by the optical flows computed for every camera. The energy to be minimized for the regularization of the set of 3D point is: Lecture Notes in Computer Science 5 E¡¯ X0, .., ¯ XNc−1¢= Nc−1 X c=0 ÃZΩ dist(¯ Xc,¯ Rc)2+α 2 X i=0 ZΩ¡∇Xc,i¢TD(¯ hc)∇Xc,i! (3) where αis a parameter that states the balance between the two terms and dist(¯ Xc,¯ Rc) denotes the distance from point ¯ Xcto the straight line ¯ Rcand is given by formula dist(¯ Xc,¯ Rc)2= 2 X i=0 ¡Xc,i −Fc,i¢2−Ã2 X i=0 lc,i ¡Xc,i −Fc,i¢!2 .(4) Our unknowns are ¡¯ X0, ..., ¯ XNc−1¢for cameras 0, ..., Nc−1 respectively. We suppose every camera has its own set of 3D points but these sets are referred to a common global set of 3D points. For consecutive cameras a large part of their set of points are going to be coincident. In particular the functions Xc are related by (2). We will see later on Sect. 4 that formulating the energy in this way allows us to create a numerical scheme by summing the contribution of every camera. The first term of (3) minimizes the distance from the 3D point ¯ Xcto the straight lines generated by the corresponding camera points passing through the camera focus. The second is a smoothness term that minimizes the variation of the 3D points according to the information given by the gradient of ¯ hc +and ¯ hc −. In this case we minimize the variation of the 3D surface by using an operator which is very similar to the Nagel–Enkelmann operator. This operator allows the method to regularize isotropically when the 3D points set varies smoothly and anisotropically when there is a strong variation of the disparity map, thus, respecting regions of different depths. Matrix D(¯ hc) is given by the following expression: D(¯ hc) = M(¯ hc +) + M(¯ hc −) ° °∇¯ hc +° ° 2+° °∇¯ hc −° ° 2(5) where ° °∇¯ h° ° 2=k∇uk2+k∇vk2 and M(¯ h) =              ∂u ∂y 2+∂v ∂y 2−³∂u ∂y ∂u ∂x +∂v ∂y ∂v ∂x ´ −³∂u ∂y ∂u ∂x +∂v ∂y ∂v ∂x ´∂u ∂x 2+∂v ∂x 2 if λ2≥° °∇¯ h° ° 2 µ1 20 01 2¶if λ2<° °∇¯ h° ° 2 .(6) This is a projection matrix on the orthogonal space to the vector field hc +/−. This matrix is similar to the Nagel-Enkelmann operator [12] but instead of using 6 Luis Alvarez, Carmelo Cuenca and Javier S´anchez λwithin the matrix we use it as a real threshold in a non-continuous function. We also consider the optical flow as the projection space and take advantage of the information in both directions. λis a threshold that states the contour value from where we obtain an anisotropic behavior. This parameter is obtained through an isotropy fraction, s, as is explained in [1]. When s→0 the method becomes anisotropical and when s→1 the method becomes isotropical. The first problem that comes out when we discretize this energy is that in every camera, the point distribution (xc(¯ X), yc(¯ X)) is given in float precision and, therefore, it does not consist of a simple pixel grid where we may discretize the equations using the surrounding neighbors. To solve this problem we use a finite element method. 4 Finite Elements Method We are going to discretize the function Xcusing the discrete set of points {Xj}j=0,..,Np−1obtained from the sequence of 2D corresponding image points (xci(Xj), yci(Xj)) , i = 0, ..., N(Xj)−1. Since (xci(Xj), yci(Xj)) are given in float precision and they are distributed in the images in an heterogeneous way, it is not convenient to use a classical finite difference scheme. So we propose a finite element scheme that takes into account the geometry of the set (xci(Xj), yci(Xj)) in every image. Following the notation introduced in Sect. 2 let us call Np(c) the number of 3Dpoints visible from camera c. For any of these points (l= 0, .., Np(c)−1), j(l) denote the associated 3D point index in a global list of 3Dpoints. We generate, using a Delaunay triangulation (see [10]), a grid of Kc ntriangles in every camera cwith (xc(¯ Xj(l)), yc(¯ Xj(l)))l=0,...,Np(c)−1points. We will call Nt(c) the number of triangles generated on camera c. In the finite element approach the Xc,i(x, y) functions are approximated as Xc,i(x, y) = Np(c)−1 X l=0 Xc,i j(l)φc l(x, y) (7) where Xc,i j(l)are the unknowns of our system and φc j(l)(x, y) is a set of basis functions centered on the nodes ¡xc(¯ Xj(l)), yc(¯ Xj(l))¢. These basis functions are defined over each triangle Kc nof vertexes (j(l), j(m), j(n)) in the following way: Dx(τ, υ)=(xc(¯ Xj(τ))−xc(¯ Xj(υ))) Dy(τ, υ) = (yc(¯ Xj(τ))−yc(¯ Xj(υ))) φc l(x, y) = (y−yc(¯ Xj(m)))Dx(n, m)−(x−xc(¯ Xj(m)))Dy(n, m) Dy(l, m)Dx(n, m)−Dx(l, m)Dy(n, m).(8) Note that the gradient of φon the triangle is constant and is given by ∇φc l(x, y) = Ã−Dy(n,m) Dy(l,m)Dx(n,m)−Dx(l,m)Dy(n,m) Dx(n,m) Dy(l,m)Dx(n,m)−Dx(l,m)Dy(n,m)!.(9) Lecture Notes in Computer Science 7 First, we note that ¯ Xc j(l)corresponds to a 3D point visible from camera cand that such a point references a 3D point, ¯ Xj,of the global 3D point list, that is, the true unknowns are the ¯ Xjpoints. Replacing Xc,i(x, y) = PNc−1 l=0 Xc,i j(l)φj(l)(x, y) in (3) we obtain: E(©¯ Xjªj=0,..,Np−1) = Nc−1 X c=0 Np(c)−1 X l=0 dist(¯ Xc j(l),¯ Rc)2 +α Nc−1 X c=0 2 X i=0 Nt(c)−1 X n=0 Np(c)−1 X l,l0=0 Xc,i j(l)Xc,i j(l0)ZKc n³∇φc j(l)´T D(uc, vc)∇φc j(l0).(10) Deriving this energy by respect to Xi jwe obtain ∂E ∂Xi j = N(¯ Xj)−1 X m=0 2³Xcm(¯ Xj),i −Fcm(¯ Xj),i −lcm(¯ Xj),i 2 X k=0 lcm(¯ Xj),k ³Xcm(¯ Xj),k −Fcm(¯ Xj),k´! +α N(¯ Xj)−1 X m=0 Nt(cm(¯ Xj))−1 X n=0 Np(cm(¯ Xj))−1 X l,l0=0,j(l0)=j 2Xc,i j(l)ZKc n (∇φc l)TD(uc, vc)∇φc l0(11) for every i= 0,1,2. In order to find the minimum of this energy, these equations are equaled to 0. A system of 3Npequations and unknowns is thus generated. 4.1 Solving the System We are going to see how to build up the system and solve it numerically. We express the system as Au =bwhere Ais a square matrix, bis the independent constant vector and uis the vector of 3Npunknowns to be ordered in the following way u= (X0 0, X0 1, ...., X0 Np−1, X1 0, ..., X1 Np−1, X2 0, ..., X2 Np−1).(12) We will concentrate on (11) to derive both matrix Aand vector b. To compute vector bwe realize that constant values generated from (11) are only affected by the first term in the following manner bj+iNp= N(¯ Xj)−1 X m=0 ÃFcm(¯ Xj),i −lcm(¯ Xj),i Ã2 X k=0 lcm(¯ Xj),kFcm(¯ Xj),k!! .(13) 8 Luis Alvarez, Carmelo Cuenca and Javier S´anchez Both terms of (11) supply information for matrix A= (aj,j0). The first term of the energy supply the following information to the diagonal part of matrix A. For every j= 0, ..., Np−1 and i= 0,1,2 we have aj+iNp,j+iNp+ = N(¯ Xj)−1 X m=0 µ1−³lcm(¯ Xj),i´2¶.(14) Note that += in this formulation means accumulation of values. Matrix A is initialized to null values. In the process of constructing this matrix, values are accumulated in different steps to facilitate the task of dividing the algorithm. Later on we will see that the diagonal part of this matrix is affected by the second term of the energy. For the non-diagonal part of Asuch a term provides the following: aj+Np,j =aj,j+Np+ = − N(¯ Xj)−1 X m=0 lcm(¯ Xj),0lcm(¯ Xj),1 aj+2Np,j =aj,j+2Np+ = − N(¯ Xj)−1 X m=0 lcm(¯ Xj),0lcm(¯ Xj),2 aj+2Np,j+Np=aj+Np,j+2Np+ = − N(¯ Xj)−1 X m=0 lcm(¯ Xj),1lcm(¯ Xj),2.(15) If we now concentrate on the second term of the energy we realize that the contribution to matrix Ais the same for the 3 coordinates i= 0,1,2. We examine every camera cand for every camera we go through every triangle Kc n. Each triangle is composed of 3 vertexes (xc(¯ Xj), yc(¯ Xj)),(xc(¯ Xj0), yc(¯ Xj0)),(xc(¯ Xj00 ), yc(¯ Xj00 )). From every triangle we compute the following contributions to the matrix: aj,j+ = ZKc n¡∇φc j¢TD(uc, vc)∇φc j aj0,j0+ = ZKc n¡∇φc j0¢TD(uc, vc)∇φc j0 aj00 ,j00 + = ZKc n¡∇φc j00 ¢TD(uc, vc)∇φc j00 aj0,j =aj,j0+ = ZKc n¡∇φc j¢TD(uc, vc)∇φc j0 aj00 ,j =aj,j00 + = ZKc n¡∇φc j¢TD(uc, vc)∇φc j00 aj00 ,j0=aj0,j00 + = ZKc n¡∇φc j0¢TD(uc, vc)∇φc j00 (16) these contributions are the same for the 3 coordinates. Lecture Notes in Computer Science 9 Fig. 2. Stereoscopic pair for Herve’s face 4.2 Numerical Scheme If we study matrix Awe will realize that most of its values are null, therefore solving the previous system will induce a lot of unnecessary calculations and, therefore, a very slow method. We have implemented an efficient numerical scheme for the resolution of the previous system. Studying the way the system matrix Ais constructed we realize that it is a sparse matrix and most of its values are situated in symmetrical positions. Also we realize that in most cases, a single point will be surrounded by a small number of triangles, inducing, for every point on the mesh, a small number of values for the second term of the energy. To take advantage of this we express the system Au =bas α(D+H)u+Cu = bwhere Cis a matrix with the contribution of the first term of the energy, Dis the diagonal matrix with the contribution of the second term of the energy and His the contribution of the non-diagonal part of the second term. To solve this system we use a fixed point equation: u= (αD +C)−1(b−αHu).(17) 5 Experimental Results In this section we present some experiments. For the first experiment we propose a single stereo pair (Herv´e’s face). The second is a bust for which we have a complete round sequence of 47 calibrated images. 5.1 Single Stereo Pair In Fig. (2) he have a stereo pair of a face. The purpose of this experiment is to show the validity of the method when there is only a single disparity map estimated for each camera. In Fig. (3) we show the profile and front of the Herve’s reconstruction obtained through the method explained in [1]. We have only reconstructed the region belonging to the face. During the process we have to construct two triangulations to be used within the finite element method. In Fig. (4) we show the Delaunay triangulations for