A variational approach to 3D cylindrical geometry reconstruction from multiple views

Electronic Letters on Computer Vision and Image Analysis 6(2):54-66, 2008 A Variational Approach to 3D Cylindrical Geometry Reconstruction from Multiple Views Luis ´ Alvarez, Carmelo Cuenca, Agust´ın Salgado and Javier S´anchez Departamento de Inform´ atica y Sistemas, Universidad de Las Palmas de G.C., 35017, Las Palmas de G.C., Spain Received 2 October 2006; accepted 22 October 2007 Abstract Inthispaperwepresentavariationaltechniqueforthe reconstructionof3Dcylindricalsurfaces. Roughly speaking by a cylindrical surface we mean a surface that can be parameterized using the projection on a cylinder in terms of two coordinates, (l, θ), representing the displacement and angle in a cylindrical coordinate system respectively. The starting point for our method is a set of different views of a cylindrical surface, as well as a precomputed disparity map estimation between pair of images. The proposed variational technique is based on an energy minimization where we balance on the one hand the regularity of the cylindrical function given by the distance of the surface points to cylinder axis, and on the other hand, the distance between the projection of the surface points on the images and the expected location following the precomputed disparity map estimation between pair of images. One interesting advantage of this approach is that we regularize the 3D surface by means of a bi-dimensional minimization problem. We show some experimental results for large stereo sequences. Key Words: Optical flow, Stereoscopic vision, Disparity map, Camera calibration, Anisotropic diffusion. 1 Introduction The problem of 3D geometry reconstruction from multiple views has received much attention during the last years. It is straight forward to recover a 3D surface from a stereoscopic pair of images since the pixel grid makes it easy to establish the relation between the points. Nevertheless when we have more than two or three images and try to recover the original scene, the problem becomes much more complex. There have been different strategies to overcome this problem. The most traditional approach is to first obtain a set of scattered points and try to find the surface that best fits the set like the one proposed by Hong-Kai et al. in [11], where authors propose a technique to find out a surface from an unorganized set of points. In [3] the authors propose a volumetric method where the 3D surface is an isosurface of the volumetric grid. Another approach are the particle-based methods like in [26] and [9] Correspondence to: [email protected] Recommended for acceptance by Jon Sporring ELCVIA ISSN:1577-5097 Published by Computer Vision Center / Universitat Aut`onoma de Barcelona, Barcelona, Spain L. ´ Alvarez et al. / Electronic Letters on Computer Vision and Image Analysis 6(2):54-66, 2008 55 in where the points have attraction and repulse forces following the Newtonian dynamics. Kanade et al. [18] proposed a deformable mesh representation to match multiple dense depth stereo data. More recent methods like [15] combines the 3D and 2D image information and in [20] the authors propose a method for uncalibrated cameras that combines different aspects from projective reconstruction, self-calibration and dense depth estimation. Some new energy minimization techniques have been proposed in the literature like the one proposed by Faugeras and Keriven in [6, 7] where a level set approach is proposed to minimize a surface energy. Other energy minimization methods are [21] and [13] which utilize graph cuts to minimize the energy functional. Normally after reconstruction the 3D surface is noisy. The previous approaches consider implicit or explicitly some regularizing means to reduce the surface noise. There are techniques that focus only on the regularization of the mesh. In papers like [24, 25, 12] the authors propose mesh smoothing techniques that in the case of the last two ones there is a mechanism of feature preserving by means of anisotropic diffusion. In this paper we propose a different approach which is also based on a variational formulation but only using disparity map estimations between pair of images. We assume that the 3D surface we want to recover has a cylindrical geometry, that is, it can be expressed as an application S: (l, θ)→R3, where (l, θ)represents a cylindrical parameterization of the 3D surface. Of course, this is an important limitation in terms of the surface geometry, but it simplifies in a strong way the complexity of the problem and it can be applied in a lot of situations like, for instance, human face reconstruction. There is an important reduction on the complexity since it is easier to estimate a unique 3D surface by projecting the 3D points in a cylindrical structure than directly from a cloud of scattered 3D points. We will also assume that the cameras are calibrated (see [5], [8] or [10] 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. The main contribution we propose in this paper is a variational model to recover the 3D geometry of a cylindrical surface. This variational model is based on the minimization of an objective function. The proposed objective function is a balance between 2 terms. In the first term we minimize the distance between the projection of the 3Dsurface points in the image sequence and the expected location following a precomputed disparity estimation between pair of images. In the second term, we regularize the cylindrical function given by the distance between the cylinder axis and the surface points. This objective function also enables regularization by preserving discontinuities on the cylindrical function. The regularizing term is similar to the terms used in other fields like stereoscopic reconstruction [1] or optical flow estimation [2]. The associated Euler-Lagrange equation of this objective function yields to a nonlinear partial differential equation that is then embedded into a gradient descend method to look for the solution. We develop an explicit numerical scheme based on finite differences to implement the method. We also present some experimental results to evaluate the method. The first experiment is composed of 36 images taken around a synthetic cylinder. For the second experiment we use a large sequence of 47 cameras located around a woman bust. The organization of the paper is as follows: In section 2 we introduce the cylindrical coordinate system necessary for the representation of the cylindrical function and the relation with the projective camera model. We also introduce the disparity map technique we use to estimate the disparity between pair of images. In section 3 we present the variational model, and we compute the associated Euler-Lagrange equations. In section 3.2 we introduce an explicit numerical scheme to discretize the model. In section 4 we present the experimental results for synthetic and real image stereo sequences and finally in section 5 the conclusions. 56 L. ´ Alvarez et al. / Electronic Letters on Computer Vision and Image Analysis 6(2):54-66, 2008 Figure 1: Cylindrical and Cartesian coordinate systems 2 The Cylinder Structure and Disparity Map Estimation 2.1 The Cylindrical Coordinate System and the Projective Camera Using the notation represented in Fig. 1 we note by N1,N2and N3the orthogonal coordinate system and by Q0the origin of the system. N1denotes the cylindrical axis. The cylindrical coordinates are expressed by means of a list of three coordinates (l, θ, r)where lis the displacement on the cylindrical axis N1,θis an angle (as it is outlined in Fig. 1) and ris the distance from a 3D point to N1. A cylindrical surface S(l, θ, r)will be given by a cylindrical function r(l, θ)in the following way : S(l, θ, r) = Q0+lN1+r(l, θ) (N2cos θ+N3sin θ).(1) S(l, θ, r)is a function that transforms a cylindrical function r(l, θ)to a surface in the Cartesian coordinate system. We will see below that our method makes use of disparity maps computed between pairs of stereoscopic images to constraint the regularization of the cylindrical function. The disparity maps are expressed in image coordinates associated to every camera. We assume the projective model for the cameras. In our problem we have Ncdifferent projective cameras and every camera is represented by a projection matrix Pcof dimensions 3x4that projects 3D points into the image plane. In projective coordinates the points in the cylindrical surface are projected onto the image plane as: ˜mc(l, θ) = Pc(S(l, θ, r),1)t.(2) To obtain the image coordinates associated to the 3D point we divide by the third component of vector ˜mc(l, θ). 2.2 Disparity Map Estimation between Pair of Images In order to estimate disparity maps between pair of images we use the technique explained in [1]. In this paper the disparity map is parameterized by a scalar function, λ, that represents the displacement on the epipolar lines. The solution for λis obtained through the minimization of the following objective function: E(λ) = ZΩ (Il(x, y)−Ir(x+u(λ(x, y)), y +v(λ(x, y))))2dx dy + C ZΩ ∇λtD(∇Il)∇λdx dy (3) where (u, v)tare the components of the optical flow between the left and right images and depend on the λ function. Ωis the image domain, Cis a positive constant, and ∇λtD(∇Il)∇λdetermines the regularization term. L. ´ Alvarez et al. / Electronic Letters on Computer Vision and Image Analysis 6(2):54-66, 2008 57 D(.)is a regularized projection matrix in the direction perpendicular to ∇Iand was first proposed by Nagel and Enkelmann in the context of optical flow estimation theory [16, 17]. The shape of this matrix is D(∇Il) = 1 |∇Il|2+ 2ν2  "∂Il ∂y −∂Il ∂x #" ∂Il ∂y −∂Il ∂x #t +ν2Id   .(4) In this formulation, Id denotes the identity matrix. When derived, we obtain a partial differential equation (PDE) that is then solved by means of a gradient descend method. To recover large disparities the method is embedded into a linear scale–space strategy. 3 The Variational Method to Recover 3D Cylindrical Surfaces The method we propose is based on the minimization of an objective function that depends on the cylindrical function. It is minimized to obtain the Euler-Lagrange equations and a gradient descend approach is then developed in order to solve the system. This gradient descend method is finally implemented by means of an explicit numerical scheme. 3.1 Minimization of the Objective Function The regularization of the cylindrical function r(l, θ)is equivalent to regularize the cylindrical surface S(l, θ). We propose a variational formulation to look for the regularized solution. This solution is the result of a minimization problem. Our model is composed of two terms: an attachment term that makes use of the disparity maps to constraint the process; and a regularizing term that is used to obtain a smooth solution. This term is designed to regularize the surface by preserving the discontinuities of the cylindrical function which are related to the varying depth of the 3D surface. The objective function proposed is E(r) = N X c=1 Z Z  mc+1 (l, θ)−mc(l, θ)−hc +(mc)  2dldθ + N X c=1 Z Z  mc(l, θ)−mc+1 (l, θ)−hc −(mc+1)  2dldθ! +αZ Z φ(k∇rk) dldθ . (5) mc+1 (l, θ)is the image coordinate for camera c+ 1 denoted by (2) and mc(l, θ)is the correspondent for camera c. Vectors hc +/−(mc) = u+/−(mc) v+/−(mc)!represent the optical flow estimations for pixel mcon camera c. Sign +corresponds to the optical flow from camera cto c+ 1 and sign −to the optical flow from camera c to camera c−1.We use a balance parameter, α. In the first term, the surface is constraint with the information supplied by the original optical flows in both senses. The second term is the regularizing term. The purpose of this term is to find a smooth solution respecting the cylindrical function discontinuities. This kind of term has been extensively used in other related subjects like optical flow estimation, stereoscopic vision or image restoration. In papers [4, 14] there is a discussion about function φ(k∇rk)and some references to other related works. A functional variation of this energy leads to the Euler–Lagrange equation that is given by the following PDE: 58 L. ´ Alvarez et al. / Electronic Letters on Computer Vision and Image Analysis 6(2):54-66, 2008 ∂E(r) ∂r = Nc X c=1 mc+1 −mc−hc +(mc)t ·∂mc+1 ∂r −∂mc ∂r − J hc + ∂mc ∂r  + Nc X c=1 mc−mc+1 −hc+1 −(mc+1)t ·∂mc ∂r −∂mc+1 ∂r − J hc+1 − ∂mc+1 ∂r  −α·div φ′(k∇rk) k∇rk∇r= 0 (6) where Jh=J u(x, y) v(x, y)!= ∂u ∂x ∂u ∂y ∂v ∂x ∂v ∂y !. In order to search for the solution we implement a gradient descend method in the way ∂r ∂t =−∂E(r) ∂r , thus having an equation in the form: ∂r ∂t =α·div φ′(k∇rk) k∇rk∇r − Nc X c=1 mc−mc+1 −hc+1 −(mc+1)t ·∂mc ∂r −∂mc+1 ∂r − J hc+1 − ∂mc+1 ∂r  + Nc X c=1 mc+1 −mc−hc +(mc)t ·∂mc+1 ∂r −∂mc ∂r − J hc + ∂mc ∂r  .(7) In this case the divergence term is well known and acts like a diffusion scheme. If we expand the divergence expression we obtain div φ′(k∇rk) k∇rk∇r=φ′(k∇rk) k∇rkrξξ +φ′′ (k∇rk)rηη (8) where η=∇r k∇rkand ξ=η⊥are the unitary vectors in the directions parallel and perpendicular to the gradient, respectively. Playing with function φ(s)it is possible to achieve an anisotropic diffusion at contours. The first in proposing this kind of diffusion equation were Perona and Malik [19] in where they introduced a decreasing function to avoid diffusion at contours. Other related workshave appeared and some of them have addressed the problem of regularizing along the isophotes. Our objective is not only to diffuse isotropically at homogeneous regions but also to diffuse along the contours. We can get this behaviour by adapting the coefficients of rξξ and rηη. The first allows the diffusion to be carried out along the isophotes and the second in the direction of the gradient. 3.2 Numerical Scheme In this section we study how to implement an explicit numerical scheme for this method. We derive ∂m ∂r analytically from (2). Regarding (8) the divergence is divided in two terms and the values for both of them are given by the following expressions: L. ´ Alvarez et al. / Electronic Letters on Computer Vision and Image Analysis 6(2):54-66, 2008 59 rξξ =rxxr2 y−2rxryrxy +ryyr2 x r2 x+r2 y , rηη =ryyr2 y+ 2rxryrxy +rxxr2 x r2 x+r2 y .(9) The first and second derivatives on xand yare approximated by finite differences as rx≃ri+1,j −ri−1,j 2, ry≃ri,j+1 −ri,j−1 2, rxx ≃ri+1,j −2ri,j +ri−1,j, ryy ≃ri,j+1 −2ri,j +ri,j−1, rxy ≃ri+1,j+1 +ri−1,j−1−ri−1,j+1 −ri+1,j−1 4.(10) The derivatives of the components of the optical flow, ∂u ∂x ,∂u ∂y ,∂v ∂x , and∂v ∂y , have also been approximated by finite differences. The final numerical scheme is implemented by means of an explicit scheme in the following way: rt+1 =rt+dt ·(α(rξξ + g (k∇rk) rηη) − Nc X c=1 mc−mc+1 −hc+1 −(mc+1)t ·∂mc ∂r −∂mc+1 ∂r − J hc+1 − ∂mc+1 ∂r  + Nc X c=1 mc+1 −mc−hc +(mc)t ·∂mc+1 ∂r −∂mc ∂r − J hc + ∂mc ∂r !   .(11) Function g(s)is a decreasing function that disables isotropic diffusion for big values of the gradient. This function depends on a parameter called λthat represents the value from where the gradient is considered to be a contour. This parameter is estimated by means of an isotropy fraction sthat states the value of the gradient from a percentage (see [1] for more details). The image coordinates, m, are computed from the cylindrical coordinates, r, θ and l. For this we have to know from what cameras the values of the cylindrical function, r(θ,l)were obtained. 3.3 Surface Initialization. Building the Initial Cylindrical Function In the proposed numerical scheme, we need to provide an initial guess for the 3D-cylindrical surface. We will use a simple technique based on the 3D surface reconstruction obtained using pair of images. Initially, we have as many 3D surfaces as cameras are there in the system. We will transform the 3D surfaces into a unique cylindrical function that will be used in the method. Given the set of 3D surfaces the main steps for computing the initial cylindrical function are: 1. Estimate Q0,N1,N2and N3. This step depends on the position and orientation of the surface on the scene. For instance if we assume that the image sequence is obtained by turning a camera around the object in a orthogonal plane to the cylinder axis, then we can easily compute the cylinder parameters as it is explained below. 60 L. ´ Alvarez et al. / Electronic Letters on Computer Vision and Image Analysis 6(2):54-66, 2008 Figure 2: Views 0, 4, 28 and 32 of the synthetic cylinder sequence. 2. Built the cylindrical image. The cylindrical function will be represented through an image. This is what we call the cylindrical image. The rows and columns of this image are given by the N1axis and the angle, θ, respectively. The bigger this image is the more accurately the 3D surface is represented. 3. Create the cylindrical function, r(θ, l). We project all 3D surfaces obtained from each camera in the cylindrical image. We merge the information of all the 3D surfaces by computing the mean for all coincident 3D points that are projected on one pixel of the cylindrical image. We also keep track of the cameras from where the 3D points are seen. The first step is to estimate the position, Q0= (q0, q1, q2)T,and axis, N1= (n11, n12, n13)T,N2= (n21, n22, n23)Tand N3= (n31, n32, n33)T,of the cylindrical coordinate system. We have supposed that the camera configuration system is cylindrical in the sense that all the cameras are situated around the scene and looking at the centre. We also suppose that the foci of the cameras are situated close to a common plane. Q0 is estimated as the average of the 3D points of all surfaces. N1is the cylindrical axis and is computed as the orthogonal direction to the plane where the camera foci lay, N2is the unitary vector in the camera focuses plane that points from the cylindrical axis to the focus of the first camera and N3is orthogonal to N1and N2. In the second step we are concerned with the problem of representing the cylindrical function through a bi-dimensional image. We have to compute the dimensions of an image that will allocate the values of the 3D points in cylindrical coordinates. To calculate the number of rows the lowest and highest 3D points in the N1component are computed. The difference between them defines the size of the cylindrical axis. The number of columns is estimated knowing that 2·π·radius is the length for the cylinder. We adapt the value of radius in order to obtain an image with regular pixels (same pixel height and width). This value depends on the dimension of the image in the N1axis. This image represents the r(θ, l)function. The last step consists of assigning a value to every pixel on the image. This process is carried out by representing the 3D points in cylindrical coordinates and computing a mean for coincident points on a pixel. There may be some locations where no 3D point is projected, so a post-processing to fill these gaps is necessary. These are filled from the values of the surrounding pixels. At the same time that we compute a value for every pixel we keep the information of the cameras from where these points were obtained. This will be necessary for the resolution of the method. Finally what we obtain is an image that for every pixel position there is a value for r(θ, l). This is the way we simplify the set of 3D surfaces into a cylindrical function that we use as an initial guest in the proposed numerical scheme. 4 Experimental Results In this section we show the results of reconstructing and regularizing a synthetic cylinder and a bust sequence. In the web page ”http://serdis.dis.ulpgc.es/˜jsanchez/research/demos” the reader can find more details on the numerical experiments we present here. In particular you can find vrml models of the 3D reconstruction. In Fig. 2 we show several views of the cylinder. This sequence is composed of 36 images around the figure. L. ´ Alvarez et al. / Electronic Letters on Computer Vision and Image Analysis 6(2):54-66, 2008 61 Figure 3: Different views of the 3D reconstruction of the cylinder. Front view, upper view, and view from the bottom. Figure 4: Regularization of the synthetic cylinder with λ= 0,1. We first compute the disparity maps between every pair of consecutive images, then we obtain a 3D surface from every disparity map and finally we build the cylindrical function. In Fig. 3 we show the final 3D cylindrical surface. In this figure we show three views: a front view, an upper view and a view from the bottom. We may appreciate that the surface is not regular. In Figs. 4 and 5 we show different regularizations. For the first example we have used a value for αof 0.1 and for the second a value of 0.5. The 3D model obtained for the second is more regular. The result is in both cases a surface with less noise. Next we show a real sequence of a bust. A lot of research has been devoted to human body reconstruction from image sequences (see for instance [22]). In this case the sequence is composed of 47 images taken around a bust. Figure 6 shows the configuration of this sequence with the projection planes of the cameras. In Fig. 7 we show different views of the sequence. The Bust is introduced in a platform with some squares and little labels on them. These labels were introduced to calibrate the cameras. From Figs. 10, 11 and 12 we may appreciate several regularizations for different values of αand s. Figure 5: Regularization of the synthetic cylinder with λ= 0,5. 62 L. ´ Alvarez et al. / Electronic Letters on Computer Vision and Image Analysis 6(2):54-66, 2008 Figure 6: Bust configuration: This figure shows the 3D reconstructed bust and the distribution of the projection planes corresponding to the 47 cameras Figure 7: Images 0, 12, 20 and 32 of the Bust sequence 5 Conclusions In this paper we have presented a novel method for the reconstruction and regularization of cylindrical surfaces. This kind of surfaces allows us to develop specific techniques which are simpler and easier to implement than other more general methods. In particular we have considerably simplified the process of 3D reconstruction from multiple views. First, we obtain a set of 3D points for each pair of stereoscopic views using a previous technique. Then, all the 3D points are merged in a common cylindrical function, which is a bi-dimensional array, in a very simple way and it is very fast since the operations involved are projections and averaging of the 3D points. This avoids the burden of determining the 3D geometry directly from the set of 3D points. We have taken advantage of the simplicity of cylindrical coordinates to represent the set of 3D points. Once the cylindrical function is built the problem of surface reconstruction and regularization is reduced to a 2D Figure 8: The left image represents the texture of the Bust sequence projected on a cylindrical image. The right image is the cylindrical function represented in grey levels (the white colour is associated to the highest values)