SciELO - Scientific Electronic Library Online

vol.13 issue4Letter from the editor-in-chiefModel-based evolution of collaborative agent-based systems author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand




Related links


Journal of the Brazilian Computer Society

Print version ISSN 0104-6500

J. Braz. Comp. Soc. vol.13 no.4 Campinas Dec. 2007 



Regularized implicit surface reconstruction from points and normals



B. MederosI; M. LageII; S. AroucaII; F. PetronettoII; L. VelhoI; T. LewinerII; H. LopesII

ILaboratório Visgraf Instituto de Matemática Pura e Aplicada {boris, lvelho }
IILaboratório Matmídia - Departamento de Matemática. Pontifícia Universidade Católica do Rio de Janeiro {sueni, mlage , fbpetro, tomlew, lopes }




We consider the problem of surface reconstruction of a geometric object from a finite set of sample points with normals. Our contribution is to present a new scheme for implicit surface reconstruction. Similarly to the multilevel partition of unity (MPU) method we hierarchically divide the domain obtaining local approximation for the object on each part, and then patch all together obtaining a global description of the object. Our new scheme uses ridge regression and weighted gradient one fitting techniques to get better stability on local approximations. The method behaves reasonably on sparse set of points and data with holes as those which comes from 3D scanning of real objects.

Keywords: Computer graphics, Implicit modeling, Surface reconstruction, Partition of unity, Ridge regression.




Surface reconstruction is playing an important role in Computer Graphics [6]. Reconstruction is a very complex problem not only because the adjacency and proximity relations of data are unknown, but also because there are a lot adversities that need to be faced. The data in which these algorithms are applied come usually from 3D scanner. Currently this devices are able to handle real objects with increasing complexity and the resulting point clouds of the data acquisition contain fine details, rapid geometric variations, complex topology and sharp features. However, the process of capturing the point cloud introduces sparse samples, holes (due to the occlusion of certain part of the object by other ones) and noise.

Several techniques have been studied to solve the surface reconstruction problem along the last two decades producing a variety of algorithms. Some methods are based on Delaunay triangulation concepts, among them, we can cite the works [13, 14, 15] for clean data and the works [12, 21] for noisy data. Other solutions are based on local parameterizations of shape like the Moving Least Square (MLS) approach [20,10,11]. Finally, there is also an important class of methods based on implicit function approximations. On this class, there are ones that use radial based functions (RBF) [8, 22, 25] and others based on domain decomposition schemes [17, 16, 19, 24, 3].

In this work we will focus on an implicit based scheme for surface reconstruction. Implicit Surface is a very useful representation for 3D objects, mainly because the inferred shape is computed by a formula which allows the computation of basic modeling operations in a relatively easy way [23]. Most of the boundary of man made objects are compose of several patches which can be approximated by algebraic surfaces. When the object's shape is complex, a common procedure is to elevate the algebraic degree in order to obtain more precision on the approximation. However, in this case due to the ill posedness of this fitting some spurious connected components appear on the reconstructed surface [5]. An alternative solution is to decompose the domain hierarchically in compact parts and obtain local approximation for the object in each part, and then patch all together in order to obtain a global description of the object. A practical scheme that uses such solution is the method of multilevel partition of unity implicit (MPU) [17]. It provides an adaptive error-controlled approximation of the signed distance function from the surface.

Problem description. Given a finite set of points in sampled from a surface S in the space, also assume that for each point p the unitary normal vector n to the surface at p is given. The objective of this work is to obtain an implicit function F : such that the isosurface of level zero, F-1(0), approximates adaptively S using local error-control.

Contributions. The reconstruction scheme proposed in this paper is an extension to 3D surface reconstruction of the work [18] proposed to tackle the 2D case of curve reconstruction. In this extension we introduce a weighted local algebraic approximation which produces a substantial improvement compared to the original two dimensional method. Our algorithm is based on the main ideas of MPU [17] that hierarchically subdivides the domain in several parts and later computes local shape approximations on each part. Our method presents a different strategy for the local shape approximations and improves its numerical stability by the use of the ridge regression technique and weighted gradient one fitting. As a consequence, it avoids the generation of spurious connected components on the reconstructed implicit surface.

Paper outline. Section 2 introduces some basic concepts. Section 3 describes the main ideas of the original MPU method. Section 4 describes some schemes for implicit surface fitting. Section 5 introduces our new method and the implementation details. Section 6 shows some results. Finally, Section 7 concludes and proposes future directions of this work.




A subset is called an Implicit Surface if there is a function , and a real number c R such that O = F-1(c). The implicit surface F-1(c) is regular if F is differentiable and satisfies the condition that at each point x F-1(c) the gradient of F at x does not vanish.

A polynomial of degree d defined on is a function Pd : given by the following expression:

An algebraic surface of degree d is the implicit surface Pd-1(0). It is convenient to adopt a suitable notation for Pd. We adopted the vectorial one proposed by Tasdizen et al. in [18], which is:




The elements of the vector a l are the coefficients ai,j,k (0 < i, 0 < j, 0 < k, and 0 < i + j + k < d) of Pd, and the elements of the vector v(x,y,z) l are the monomials of Pd. The dimension l of the vectors a and v(x,y,z) depends on the degree d and is obtained by the following expression:


A partition of unity (PU) [1],[4] is a mathematical tool very useful to combine local approximation in order to define a global one. Important properties such as the global maximal error and the convergence order could be inherited from the local approximation. The basic ideas behind the global approximation construction using partition of unity are the following:

i.  divide the domain in parts,

ii.  obtain a local approximation for each part using a subset of the data that belongs to it,

iii.  obtain a global approximation by the use of a weighted combination of local solution through the use of smooth non-negative functions that correspond to the weights. In each point of the domain, the sum of these weight functions should be one.

More precisely, consider a compact domain Ω and denote by {φi}i=1,...,n the set of non-negative functions with compact support such that:

Let be a set of functions defined in supp(φi). Each function in represents a local approximation for the points of that belong to supp(φi). A global approximation for the function f : Ω could be obtained as following:

where fi . Consider {wj}j=1,...,n a set of non-negative functions with compact support such that:

The partition of unity functions φi could be generated by the following equation:

The main idea of partition of unity could be resumed by the equations (3) and (4). Such equations form the basis of the algorithm Multilevel Partition of Unity Implicit (MPU) proposed by Ohtake et al. in [17].



The method called Multilevel Partition of Unity (MPU) was proposed by Ohtake et al. in [17] originally to build an implicit surface approximation of a set of points and normals in . The MPU uses a partition of unity to obtain a global implicit surface approximation for the boundary of the object combining local approximations. It uses an octree as an hierarchical scheme to guide the domain subdivision.

Follows a concise description of how the MPU builds an implicit function that globally approximates the points.

The method initially centers the point of at the origin. After that, the points are scaled in such a way that the square Ξ = [-1,1]3 contains all points of . We will adopt the same name for the set of points after these two transformations.

The method builds an octree by the use of recursive procedure where the subdivision of each node is controlled by the error of the local approximation. In other words, the refinement criteria for a node i of the octree consists of computing the local error of the approximation and when this error is greater than a given tolerance, then the node is subdivided in eight new nodes and recursively the same test is again used for each one of its child nodes.

Each node i on the octree is associated to a weight function wi with a compact support that is used for the partition of unity global approximation. The compact support of wi is defined as a circle of radius ri centered on the middle of the node i. Such radius is chosen proportionally to the size of the diagonal of the square corresponding to the node i, denoted by di.

In the MPU original method, a quadratic spline function b : is used to build the weight function wi:

where ci is the center of the corresponding node i in the octree. The value of wi is zero outside the support region. According to the equation (4), the partition of unity function φi associated to the node i is defined as:

where nl is the number of leaves on the octree.

The MPU method uses a quadric (degree 2 polynomial function) to locally approximate the signed distance function to the boundary of the object. According to the distributions of the normal in a node of the octree either a 3D quadric or bivariate quadric polynomial are used to approximate the local shape. To find the coefficients of the quadric function Qi : supp(φi) for the node i, Ohtake et al. [17] used a least squares scheme that will be reviewed in subsection 4.1.

At each node, a least squares problem is solved considering only the points on that belong to its support region as input. Sometimes (especially when the density of is not uniform) the circle of radius ri of a node i doesn't contain a sufficient number of points to estimate robustly the quadric that approximates such points. If the number of points on the support region of a node is not sufficient to solve the minimization problem, then they adopted a solution that increases the radius of the support region until such minimal condition is guaranteed.

They suppose that the surface S from where the points are sampled is a level 0 isosurface of a function f : . Then, they can use the partition of unity equation (7) to obtain a F : that globally approximates f:

For more information of this method, see [17].



Follows the description of some methods to obtain an implicit surface approximation.


Suppose that a set of q points = {p1,p2,...,pq} sampled from a 3D surface S is given. A simple way to obtain an algebraic surface Pd-1(0) that approximates S is to minimize the total sum of the squared algebraic distance, denoted by ealg, from each point pi to the curve Pd-1(0). Using the vectorial representation of Pd in (1), we can write ealg as:

where vi corresponds to the evaluation of the vector v, see equation (2), at the points pi = (xi,yi,zi). This minimization is subject to the constraint that the sum of the squared coefficients of Pd is equal to one (||a||2 = 1). Note that without this constraints the minimum of ealg is reached at the null vector. A more robust approximation measure is to consider the weighted algebraic distance

where wi = w(pi) > 0 is a set of weight and w(x, y, z) is a positive weight function such that w(x) approaches to zero when ||x|| goes to infinity or a compactly support function. To improve the notation, define the matrices M of size l × q and S of size l × l as follows:


Thus, the minimization problem we have to solve is the following:

By the use of the Lagrange multiplier λ, the constrained minimization problem is rewritten as:

Its solution is the unitary eigenvector of S associated to the eigenvalue of S with smallest value [7].

Although this method is invariant to affine transformations [7], it has some delicate problems. Its results are sensitive to small perturbations on the input data. Moreover, the algebraic surface Pd-1(0) doesn't consider the points continuity. Thus, it could generate undesired connected components or glue components that are originally separated. For more details see [2, 18, 7].


To avoid the problems of continuity and of sensitivity caused by small perturbation on the input data [2] proposed a new method that considers not only the set P but also a given set of normal unitary vectors N = {n1,n2,...,nq}, where the vector ni is the normal vector to the surface S at the sampled point pi .

The gradient

of Pd at the point pi = (xi, yi, zi), will be used by this method to approximate the direction of the given normal. It is important to notice that if the gradient of Pd at pi = (xi, yi,zi) is not equal to zero, then it is perpendicular to the tangent plane of the level surface Pd that pass through the point pi.

The gradient one fitting (GOF) method [2] is in fact an unconstrained least square problem. We have introduced a modification in the original GOF method that consist of a weighted least square problem as follows

where μ is the weight given to the terms that approximates the normals and wi is a set of weight depending of each pi. In order to follow the vectorial representation of Pd, we define the following matrices and vectors:

• The matrix of size l × 3:

• The gradient vector Pd:

• The matrix SN of size l × l:

• The vector gN of size l:


Therefore, the weighted GOF optimization problem can be rewritten as:

And its solution is obtained by solving the following system of linear equations:


The MPU 3D fitting is based on the minimization of the squared algebraic distance eaíg (8) plus an energy term empu which in some way tries to favor a local approximation of the signed distance function to by 3D quadric.

The empu depends on the set of auxiliar points Q = {q1,q2,...,qm} which is a subset of the vertices of the node in the subdivision octree, and this implies that the maximum cardinality of Q is eight. Following [17] a vertice q of the node is reliable for computing the approximated average signed distance, (where p1, p2 and p3 are its three near sample points on ) if n1(q - p1), n2(q - p2)and n3(q - p3) have the same sign.

Then the energy empu is defined as

where dj is the approximated average signed distance of qj to . The combination of empu with GOF is a good way to determine local shape approximation which takes advantage of the benefits of (9) and (11)

The minimization of (12) conduces to the following linear optimization problem

where SMPU and gMPU are defined as follows:


The minimization problem (13) conduces to the following linear system


When the matrix doesn't have a maximal rank or is ill conditioned then the technique called ridge regression (RR) can be used. Statisticians use it frequently to remove the collinearity of the input data. The first proposal to obtain algebraic surfaces that fits better was done by Tasdizen et al. in [18]. The RR technique basically modifies the optimization problem of the gradient one fitting method by adding a new term:

where Δ is a diagonal matrix of size l × l and the real constant κ determines the weight given to the new term. The minimization problem solution is obtained by solving the following system of linear equations:



We propose in this section a new method for implicit surface fitting. It combines the MPU scheme and the gradient fitting one and ridge regression methods to improve the implicit surface approximation.

Data input. Our algorithm considers as input data a set of q points , and a set of the corresponding q unitary normal vectors N.

We are assuming that the points of have been translated in such a way that the center of mass is the origin of the coordinate system, and also that they have been scaled in such a way that all points are contained in the square Ξ = [-1,1]3. Such square Ξ is the starting region for the hierarchical adaptive space subdivision guided by the use of an octree data structure.

Support regions. We use the same octree support region scheme of the MPU method. The support region for each node i is the disk of radius ri = adi centered at ci.

Local approximations. In our method, we adopt the mpu 3d fitting, gradient one fitting and ridge regression techniques to obtain the coefficients of a degree d algebraic function Pdi for the local approximation at the node i. The local approximation is only computed when the node contains sample points on its support region. To do this task, we consider the set of points i = {pj1,...,pjqi} that are on the support region of the node i, in case that we have enough points (a number bigger than Nmin = 15), otherwise we grow the support region until we get a minimum number of points. To run the RR method, we have also to equip each point of i with its unitary normal vectors. Thus, using the same indexes of these vectors we construct the sets Ni N for each node i.

To determine the local approximation to the shape of the surface on the node i we take a similar strategy to the MPU algorithm, i.e, two types of local approximations are computed: 1) a 3D quadric, 2) a bivariate quadratic polynomial in local coordinates. The first is used to approximate parts of the surface composed of more than one surface sheet and the second to approximate a local smooth patch. To determine which of these kinds of approximation we will compute we apply the same method proposed in [17] which roughly speaking is based on computing an average normal direction on Ni, if the maximum deviation of the normals to the average normal nave is bigger than π/2 then we compute 1) otherwise we compute 2), for more detail, see [17].

In the first case the coefficient of 3D quadric are determined minimizing the objective function:

where the matrices Si, SN,i, SMPU,i and the vectors gN,i, gMPU,i are computed for the node i using the expressions presented in Section 4. However, we consider as input for this computation the sets Pi and Ni. The solution to the minimization problem is obtained by solving the following system of linear equations:

In the second case a bivariate quadratic polynomial is determined with the domain being the plane Πi orthogonal to nave passing through the centroid of the points in i minimizing


and are the coordinates of the points pj i in an orthogonal system of coordinates (u, v, nave) where u,v is an orthogonal base in the plane is a bivariate polynomial of degree d, the minimization of (17) conduces to solving a 2 × 2 linear system of the same type of (10).

Octree construction. The octree is built using a recursive procedure, whose refinement criteria is the local approximation error. Consider a tolerance e for the local error. The condition that determines whether node i of the octree at level li should be refined is the following boolean expression:

where the local approximation error ei at node i is the widely known Taubin error metric and also the mean squared algebraic distances from the points i to the obtained surface has been considered. It seems that the Taubin error metric produces good result in less time.

Parameters of the method. In conclusion, the parameters of the method are the following:

d : degree of the algebraic surface.

lmax : maximum level for the octree

α : constant that multiplies the diagonal size of the node to obtain the radius of the support disk.

μ : weight given to the GOF term on the objective function.

μ1 : weight given to the MPU term on the objective function.

κ : weight given to the ridge regression additional term on the objective function.

• є : threshold value for the refinement condition controlled by the local approximation error.

Notice that with this set of parameters we can unify several methods presented in this paper. For example: if we want to run the original GOF method we have to assign lmax = 0 and κ= 0; if we want the original RR method we have only to assign lmax = 0.

Global approximation and function evaluation. We have supposed that the surface can be written as S = f-1(0) for some function f : . Thus, by the use of the partition of unity equation (7) we can obtain a function F : that globally approximates f:



Our method has been tested with several sparse data sets, Figures 1, 3 and 4 show the reconstruction of the points sets sampled from the bunny, knot and armadillo surface respectively using our method. In all these examples the ridge regression parameter have been fixed to κ= 0.0. The surfaces were tiled using an implementation of Topological Marching Cubes [9]. The time involved is very similar to the original MPU [17]: within 1% for the optimization, and identical for the evaluation. The typical parameters of our method are μ1 = 1.0,μ = 0.01,κ = 0.001,lmax [20,30], є [0.0005,0.005] and α = 0.75.



Figure 2






In Figure 2 in the left we have the result of our implementation of the MPU algorithm on the dragon data set in which we can see a small connected component and some artifacts on the model surface, the right picture shows the results of our method (MPU + gradient one fitting + ridge regression) we gain more numerical stabillity and we are able to remove the spurious small connected component and some of these artifacts. We tested the algorithm on this data set for different values of κ ranges from 0.005,0.001, 0.0005, 0.0002 and for these values we were able to remove this small connected component and some of these artifacts, for κ= 0.0001 a small spurious component appear on the reconstructed surface. The use of ridge regression is necessary to remove the artifacts and the connected component a combination of (MPU + gradient one fitting) is not enough to get a good result.

The effect of gradient one fitting is illustrated on the squirrel model (Figure 5), with incomplete point cloud, similarly to the usual output of laser scanner: the whole bottom part and details of the eyes and on the top are missing. Subfigures c and e show the results of the only MPU method. Observe that, due to the big hole on the head, the MPU alone generated a bump and also inaccuracy on the right eye. The combination of the MPU with gradient one fitting (Subfigures d and f) is able to reconstruct the model without bump on the head and has a greater accuracy on the eyes. In this example the use of the ridge regression does not significantly alter the results with gradient one fitting. In general, the use of the gradient one fitting can be useful in the presence of small holes since it uses the neighboring normals.



We proposed a new method that combines two powerful techniques: the weighted gradient one fitting + ridge regression and the multilevel partition of unity. On one side, the ridge regression method has been considered by the pattern analysis community as one that gives a better fitting, since it tries to have a correct topology on the surface reconstruction. However, when the surface has a complex shape it is necessary to elevate the degree of the algebraic surface to get a good result. On the other side, the multilevel partition of unity is an implicit method that is now one of the most important reconstruction techniques. In order to compute local approximations, it uses a complicate objective function. Thus, our surface reconstruction scheme not only takes the advantage of these two well recognized methods, but also unifies those methods in a simple setting.

We plan to continue this work in several directions: one possible direction is to determine an approximated tangent plane T to the samples on the node of the tree and over this plane we consider the surface as a height field and determine a bivariate local approximation Q(x,y): T . Also over T we can consider better approximations like for example using wavelets method in order to be able to faithfully reproduce the oscillations and details (texture) on each region.



This work was partly supported by the Brazilian Ministry of Science and Technology (MCT/CNPq Universal 02/2006), and the State of Rio de Janeiro (FAPERJ Primeiros projetos 2004).



[1] I. Babuska and J. Melenk. The partition of unity method. International Journal of Numerical Methods in Engineering, 40:727-758, 1997.         [ Links ]

[2] M. M. Blane, Z. Lei, H. Çivi, and D. Cooper. The 3L algorithm for fitting implicit polynomial curves and surfaces to data. Transactions on Pattern Analalysis and Machine Intelligence, 22(3):298-313, 2000.         [ Links ]

[3] Y.-L. Chen and S.-H. Lai. A partition-of-unity based algorithm for implicit surface reconstruction using belief propagation. In Shape Modeling and Applications, pages 147-155. IEEE, 2007.         [ Links ]

[4] R. Franke and G. Nielson. Smooth interpolation of large sets of scattered data. International Journal of Numerical Methods in Engineering, 15:1691-1704, 1980.         [ Links ]

[5] J. P. Gois, V. Polizelli, T. Etiene, E. Tejada, A. Castelo, T. Ertl, and L. G. Nonato. Robust and adaptive surface reconstruction using partition of unity implicits. In Sibgrapi, pages 95-104. IEEE, 2007.         [ Links ]

[6] M. Levoy, K. Pulli, B. Curless, S. Rusinkiewicz, D. Koller, L. Pereira, M. Ginzton, S. Anderson, J. Davis, J. Ginsberg, J. Shade, and D. Fulk. The digital Michelangelo project: 3D scanning of large statues. In Siggraph, pages 131-144. ACM, 2000.         [ Links ]

[7] G. Taubin. Estimation of planar curves, surfaces, and nonplanar space curves defined by implicit equations with applications to edge and range image segmentation. Pattern Analysis and Machine Intelligence, 13(11):1115-1138, 1991.         [ Links ]

[8] J. Carr, R. Beatson, J. Cherrie, T. J. Mitchell, W. R. Fright, B. C. McCallum, and T. R. Evans. Reconstruction and representation of 3D objects with radial basis functions. In Siggraph, pages 67-76. ACM, 2001.         [ Links ]

[9] T. Lewiner, H. Lopes, A. W. Vieira, and G. Tavares. Efficient implementation of Marching Cubes' cases with topological guarantees. Journal of Graphics Tools, 8(2):1-15,2003.         [ Links ]

[10] M. Alexa, J. Behr, D. Cohen-Or, S. Fleishman, D. Levin, and C. Silva. Point set surfaces. In Visualization, pages 21-28. IEEE, 2001.         [ Links ]

[11] M. Alexa, J. Behr, D. Cohen-Or, S. Fleishman, D. Levin, and C. Silva. Computing and rendering point set surfaces. Transactions on Visualization and Computer Graphics, 9(1):3-15, 2003.         [ Links ]

[12] T. K. Dey and S. Goswami. Provable surface reconstruction from noisy samples. Computational Geometry: Theory and Applications, 35(1-2):124-141, 2006.         [ Links ]

[13] N. Amenta and M. Bern. Surface reconstruction by voronoi filtering. Discrete and Computational Geometry, 22(4):481-504, 1999.         [ Links ]

[14] N. Amenta, S. Choi, T. Dey, and N. Leekha. A simple algorithm for homeomorphic surface reconstruction. In Symposium on Computational Geometry, pages 213-222. ACM, 2000.         [ Links ]

[15] N. Amenta, S. Choi, and R. Kolluri. The Power Crust, unions of balls, and the medial axis transform. Computational Geometry: Theory and Applications, 19(2-3):127-153, 2001.         [ Links ]

[16] Y. Ohtake, A. Belyaev, and M. Alexa. Sparse low-degree implicit surfaces with applications to high quality rendering, feature extraction, and smoothing. In Symposium on Geometry processing, page 149. Eurographics, 2005.         [ Links ]

[17] Y. Ohtake, A. Belyaev, M. Alexa, G. Turk, and H.-P. Seidel. Multi-level partition of unity implicits. In Siggraph, volume 22, pages 463-470. ACM, 2003.         [ Links ]

[18] T. Tasdizen, J. P. Tarel, and D. B. Cooper. Algebraic curves that work better. In Computer Vision and Pattern Recognition, volume II, pages 35-41. IEEE, 1999.         [ Links ]

[19] M. Lage, F. Petronetto, A. Paiva, H. Lopes, T. Lewiner, and G. Tavares. Vector field reconstruction from sparse samples. In Sibgrapi, pages 297304. IEEE, 2006.         [ Links ]

[20] D. Levin. Mesh-independent surface interpolation, pages 37-49. Springer, 2003.         [ Links ]

[21] B. Mederos, N. Amenta, L. Velho, and L. de Figueiredo. Surface reconstruction for noisy point clouds. In Symposium on Geometry Processing, pages 53-62, 2005.         [ Links ]

[22] V. V. Savchenko, A. Pasko, O. Okunev, and T. Kuni. Function representation of solid reconstructed from scattered surface points and contours. Computer Graphics Forum, 14:181-188, 2002.         [ Links ]

[23] A. Sharf, T. Lewiner, G. Shklarski, S. Toledo, and D. Cohen-Or. Interactive topology-aware surface reconstruction. In Siggraph, pages 43.1-43.9, 2007.         [ Links ]

[24] I. Tobor, P. Reuter, and C. Schlick. Reconstructing multi-scale variational partition of unity implicit surfaces with attributes. Graphical Models, 68(1):25-41, 2006.         [ Links ]

[25] G. Turk and J. O'Brien. Modeling with implicit surfaces that interpolates. Transaction on Graphics, 21:855-873,2002.         [ Links ]

Creative Commons License All the contents of this journal, except where otherwise noted, is licensed under a Creative Commons Attribution License