US Patent No. US10417822B2
Non-Convex Hull Surfaces
by Gabriel Taubin
Granted on 09/17/2019
Download:US10417822B2.pdf
US Patent No. US10445400B2
Non-Convex Hull Surfaces
by Gabriel Taubin
Granted on 09/15/2019
Download:US10445400B2.pdf
Fast Non-Convex Hull Computation
by Julián Bayardo Spadafora, Francisco Gómez-Fernandez, and Gabriel Taubin
Proceeding of the 2019 International Conference on 3D Vision (3DV)
doi:10.1109/3DV.2019.00087
Download:Paper, Slides
3D surface reconstruction usually begins with a point cloud and aims to build a representation of the object pro- ducing that point cloud. There are several algorithms to solve this problem, each with different priors over the point cloud, such as the type of object represented, or the method by which it was obtained. In this work, we focus on an algorithm called Non-Convex Hull (NCH), which recon- structs surfaces through a concept similar to the Medial Axis Transform. A new algorithm called Shrinking Planes is proposed to compute the NCH, based on the Shrinking Ball method with a few improvements. We prove that the new method can approximate surfaces to arbitrarily small er- ror, and evaluate its performance on the surface reconstruc- tion task. The new method maintains the same reconstruc- tion quality as the Na ̈ıve Non-Convex Hull method, while achieving a large performance improvement.
Non-Convex Hull Surfaces
by Gabriel Taubin
Proceeding of SIGGRAPH Asia 2013 Technical Briefs, Article No. 2
doi:10.1145/2542355.2542358
Download:Paper,Slides,Slides (4/page)
We present a new algorithm to reconstruct approximating watertightsurfaces from finite oriented point clouds. The Convex Hull (CH) ofan arbitrary set of points, constructed as the intersection of all thesupporting linear half spaces, is a piecewise linear watertightsurface, but usually a poor approximation of the sampled surface. Weintroduce the Non-Convex Hull (NCH) of an oriented point cloud as theintersection of complementary supporting spherical half spaces; oneper point. The boundary surface of this set is a piecewise quadraticinterpolating surface, which can also be described as the zero levelset of the NCH Signed Distance function. We evaluate the NCH SignedDistance function on the vertices of a volumetric mesh, regular oradaptive, and generate an approximating polygonal mesh for the NCHSurface using an isosurface algorithm. Despite its simplicity, thissimple algorithm produces high quality polygon meshes competitive withthose generated by state-of-the-art algorithms. The relation to theMedial Axis Transform is described.
Introduction
Figure 1 A: A 2D oriented point cloud. B: A supportinglinear half space for one of the oriented points. C: The orientedconvex hull (OCH) of the point cloud. D: An outside supporting circlefor one of the points. E: Inside supporting circles are obtained byinverting the orientation vectors. F: The non-convex hull (NCH) ofthe oriented point cloud is the intersection of the complement of allthe outside supporting circles. G: A 3D oriented point cloud. NCHsurface reconstructions based on octrees of depth 7 (H), 8 (I), and 9(J).
In this paper we are concerned with the problem of reconstructing anoriented watertight surface approximating a finite set of points withassociated unit length orientation vectors, consistently oriented withrespect to the sampled surface \(S\). Oriented point clouds are producedby laser scanners, structured lighting systems, multi-view stereoalgorithms, and simulation algorithms.
The prior art on surface reconstruction from point clouds isextensive, spanning more than two decades. Most combinatorialalgorithms produce interpolating polygon meshes, and some come withguaranteed reconstruction quality [Bernardini et al. 1999; Amenta etal. 2001; Dey 2007].
Since implicit surfaces are watertight, most approximating surfacereconstruction methods produce implicit surfaces, and throughvariational formulations reduce the problem to the solution of largesparse linear systems [Kazhdan et al. 2006; Alliez et al. 2007; Man-son et al. 2008; Hoppe et al. 1992; Boissonnat and Cazals 2002;Calakli and Taubin 2011; Alexa et al. 2003; Fleishman et al. 2005;Ohtake et al. 2005]. The estimated implicit function is oftenevaluated on a regular grid of sufficient resolution, and a polygonmesh approximation is generated using an isosurface algorithm such asMarching Cubes [Lorensen and Cline 1987]. Since the cost of estimatingor evaluating the implicit function on a regular grid is oftenexcessive, some methods perform the computations on adaptivevolumetric meshes such as octrees which require more complexcontouring algorithms [Ohtake et al. 2005; Kazhdan et al. 2006; Mansonet al. 2008; Calakli and Taubin 2011].
The method proposed in this paper falls somewhere in between thesecategories. A continuous interpolating piecewise quadratic NCH SignedDistance function is constructed as a function of the oriented pointlocations and orientation vectors through a simple and directcomputation. Then the NCH Signed Distance function is evaluated on thevertices of a volumetric mesh, such as a regular voxel grid or octreeconstructed as a function of the point locations. Finally anisosurface algorithm is also used to generate an approximatingpolygonal mesh. When the volumetric mesh is conforming, the polygonmesh is guaranteed to be watertight.
Figure 2 A: An oriented point cloud with approximately 25Kpoints. B: The polygon mesh extracted by the naive algorithm on a5003 voxel grid. C: The oriented points superimposed with the mesh.D: Detail view of the point cloud showing points and orientationvectors. E: Close-up view of B. Close-up view of C. Note that thepolygon density is higher than the point cloud sampling rate. Thereconstructed polygon mesh has 556,668 vertices and 555,386 faces.
The Naïve Algorithm
To emphasize the simplicity of the proposed method, in this section wedescribe what we call the Naïve NCH Surface Reconstructionalgorithm. In subsequent sections explain why it works, andvariations. The Naïve NCH Surface Reconstruction algorithm for afinite set of oriented points comprises three steps: 1) estimating oneNCH Signed Distance parameter for each data point; 2) evaluating the NCH SignedDistance function on the vertices of a volumetric mesh such as aregular voxel grid or octree; 3) approximating the zero level setof the NCH Signed Distance function by a polygon mesh using anisosurface algorithm such as Marching Cubes [Lorensen and Cline 1987].
For a finite set of points \({\cal P}=\{ p_1,\ldots,p_N\}\), withassociated unit length orientation vectors \(n_1,\ldots,n_N\) we definethe value of the NCH Signed Distance function at a 3D point \(x\) as themaximum over \(N\) basis functions,
\begin{equation}f(x) = \max_{1\leq i\leq N} f_i(x)\label{eq:nch-signed-distance-function-finite}\end{equation}
where for each oriented point \((p_i,n_i)\),we have one basis function
\begin{equation}f_i(x) = n_i^t(x-p_i)-\rho_i\|x-p_i\|^2\label{eq:nch-signed-distance-basis-function-finite}\end{equation}
The parameter \(\rho_i\) is set equal to zeroif the set \(J_i\) of indices \(j=1,\ldots,N\) such that\(n_i^t(p_j-p_i)>0\) is empty, and otherwise
\begin{equation}\rho_i = \min \left\{\frac{n_i^t(p_j-p_i)}{\|p_j-p_i\|^2} : j\in J_i \right\} > 0\;.\label{eq:nch-rho-finite}\end{equation}
Figure 2 shows one result obtainedwith exactly this algorithm.
The main advantage of the Naïve NCH Surface Reconstruction algorithmis its simplicity, since it can be implemented literally with only afew lines of code. The main disadvantage of the method is that itscomplexity is quadratic in the number of points. We address this issuethrough an adaptive subsampling approach which yields NCH Surfaceswhich interpolate only a subset of the input points, and approximatesthe remaining points under user-specified maximum error.
Even though large areas of missing data points and holes are filledbecause the output mesh is watertight (except for its intersectionwith boundaries of the bounding box), the algorithm not always fillsholes in an intuitive manner, as can be observed in figure 1.Some of the surface reconstruction algorithms based on variationalprinciples mentioned in the introduction tend to fill holes in a moreintuitive and predictable fashion. However, since it is not iterative,this algorithm is robust, and in many cases it can deal gracefullywith uneven sampling, as shown in figure 5 .
Non-Convex Hull
In this paper we refer to a half space as a set \(H = \{ x : f(x)\leq 0\}\) of points satisfying an inequality constraint for a continuousreal-valued function \(f(x)\) defined for every point \(x\) in a certaindomain \(U\) contained in the ambient space (2D or 3D here). A linearhalf space is defined by a linear function \(f(x)\). Given a set ofpoints \({\cal P}\), finite or infinite, and not necessarily oriented,the half space \(H\) defined above (with \(f(x)\) linear or non-linear) issaid to be a supporting half space for \({\cal P}\) if the following twoconditions are satisfied: 1) the set \({\cal P}\) is contained in \(H\),i.e., \({\cal P}\subseteq H\); and 2) there is at least one point \(p\) in\({\cal P}\) where the function attains the value zero \(f(p)=0\). TheConvex Hull \(\hbox{CH}({\cal P})\) of the set \({\cal P}\)can be defined as the intersection of all the supporting linear halfspaces for \({\cal P}\). Since a linear half space is a convex set, andconvexity is preserved by intersection, \(\hbox{CH}({\cal P})\) is alsoa convex set.
For the finite set of points \(\{ p_1,\ldots,p_N\}\), with associatedunit length orientation vectors \(n_1,\ldots,n_N\) we adopt a morerestricted definition which takes into account the point orientations.Each oriented point \(p_i\) defines a linear half space \(H_i = \{ x :f_i(x)\leq 0\}\) for the linear function \(f_i(x)=n_i^t(x-p_i)\). Wedefine the Oriented Convex Hull of the set of points \(\hbox{OCH}({\cal P})\) as the intersection of all the supporting linear half spaces\(H_i\). Note that if the orientation of the points is reversed, acompletely different result is usually obtained, since the family ofsupporting linear half spaces is different, and often even empty.
Figure 3 The geometry of the spherical support functions \(f_p(x)\).
Since the Oriented Convex Hull is a convex set, it cannot approximatethe boundary surface of an object with concavities. To be able toapproximate these surfaces we need to augment the family of supportinghalf spaces. It is necessary for this family to include non-convexhalf spaces, so that their intersection can represent solid objectswith concavities. For each point \(p_i\) in the data set \({\cal P}\) withassociated orientation vector \(n_i\), and every positive value of\(r>0\), we consider the function
\begin{equation}f_i^r(x)=\frac{1}{2r}\left\{r^2-\|x-(p_i+r\,n_i)\|^2\right\}\label{eq:spherical-support-function-r}\end{equation}
defined for \(x\) in the domain \(U\). This function is positive inside asphere of radius \(r\) centered at the point \(q=p_i+r\,n_i\), negativeoutside the sphere, attains its maximum value \(r/2\) at the centerpoint \(q\), and has unit slope \(\|\nabla\!f_i^r(x)\|=1\) at every point\(x\) where\(f_i^r(x)\) is equal to zero.
As shown in figure 3, the point\(q=p_i+r\,n_i\) is located on the ray defined by the point \(p_i\) andvector \(n_i\) at distance \(r\) from \(p_i\). In addition, \(f_i^r(p)=0\),and \(\nabla\!f_i^r(p)=n_i\) for all values of \(r\). Now we define \(r_i\)as the largest value of \(r\) so that \(f_i^r(p_j)\leq 0\) for every otherpoint \(p_j\in{\cal P}\). For finite sets of oriented points we have\(r_i>0\) for each data point \(p_i\in {\cal P}\). As a result, the halfspace \(H_i\) defined by the function \(f_i(x)=f_i^{r_i}(x)\) of equation\ref{eq:nch-signed-distance-basis-function-finite} is supporting,where \(\rho_i=1/(2r_i)>0\). We refer to these sets \(H_i\) ascomplementary spherical supporting half spaces. For the linearsupporting half spaces to be included as special cases, we need toallow \(\rho_i=0\), or \(r_i=\infty\). Note that, as opposed to theOriented Convex Hull case, here every oriented point \(p_i\) in the dataset has an associated supporting half space \(H_i\), and if theorientations are reversed (\(n_i\mapsto -n_i\)), completely differenthalf spaces are obtained.
The value of \(\rho_i\) for an oriented point \(p_i\) iscomputed as the minimum over all the positive values
\begin{equation}\rho_{ij} =\frac{n_i^t(p_j-p_i)}{\|p_j-p_i\|^2}\label{eq:rho-ij}\end{equation}
over all \(j\neq i\). Note that if \(\rho_{ij}\leq 0\) for all \(j\neq i\),then we should set \(\rho_i=0\), because in this case the linear halfspace \(H_i=\{x:n_i^t(x-p_i)\}\) is supporting for the set \({\cal P}\).
We define the Non-Convex Hull of the oriented point set \({\cal P}\),denoted \(\hbox{NCH}({\cal P})\), as the intersection of all thecomplementary spherical supporting half spaces \(H_i\), as definedabove. This is the same as saying that the complement of\(\hbox{NCH}({\cal P})\) is a union of balls. Note that\(\hbox{NCH}({\cal P})\) is also a half space. Namely, the half spacedefined by the continuous signed distance function \(f(x)\) shown inequation \ref{eq:nch-signed-distance-function-finite}. This functionis well defined when the data set \({\cal P}\) is bounded, and inparticular when it is finite.
Figure 4Results on evenly sampled low noise surfaces. Left: Oriented points.Center: Reconstruction with an octree of depth 9.Right: Reconstruction with an octree of depth 10.
Relation to the Medial Axis Transform
In this section we assume that the set of points \({\cal P}\) is asampling of the boundary surface \(S\) of a bounded solid object \(O\),and that the object \(O\) is an open set in 3D. As a result, thesurface \(S\) is bounded, orientable (separates the inside from theoutside of \(O\)), closed, and it has no boundary (i.e., noholes). Furthermore we assume that it is smooth, has a continuous unitlength normal field pointing towards the inside of \(O\), and hascontinuous curvatures.
Figure 5Results on unevenly sampled surfaces. Left: Oriented points.Center: Reconstruction with an octree of depth 9.Right: Reconstruction with an octree of depth 10.
The Medial Axis Transform (MAT) is a representation of the object \(O\)as a union of balls. A ball \(B=B(q,r)=\{x:\|x-q\|
We define the Medial Axis \(\hbox{MA}(O)\) of \(O\) as the set of centersof medial balls. Since two different medial balls cannot have the samecenter, the mapping \(\hbox{MA}(O)\rightarrow \hbox{MAT}(O)\) whichassigns each medial ball center to the corresponding medial ball iswell defined, 1-1 and onto. Another way of describing the Medial AxisTransform is as a set of points called Medial Axis Set, augmented witha non-negative radius function which assigns to each medial ballcenter the corresponding medial ball radius. These radii are, ofcourse, not independent of each other.
In this paper we present an alternative description of the Medial AxisTransform, where each medial ball is not described by its center andradius, but by one of its boundary points and the radius. If \(B\) is amedial ball of center \(q\) and radius \(r\), \(p\) is a point in theintersection of the boundary of \(B\) and \(S\), and \(n_p\) is the surfacenormal to \(S\) at \(p\) pointing towards the interior of \(O\), since theboundary of \(B\) and \(S\) are tangent, the ball center \(q\) must lie onthe normal ray defined by \(p\) and \(n\), in which case we have\(q=p+rn_p\). In addition, because of the maximality of the ball \(B\),the radius \(r\) is uniquely determined: it must be equal to the maximumof radii \(r'>0\) of balls centered at points \(q'=p+r'n\) lying on theray defined by the point \(p\) and the vector \(n\), fully contained in\(O\). On the other hand, if \(p\) is a point on the surface \(S\), since\(O\) is the union of all the medial balls, and \(S\) is the boundary ofthis set, a medial ball \(B\) must exists so that \(p\) belongs to theintersection of the boundary of \(B\) and \(S\).
In summary, every medial ball can be written as \(B(p+r_p n_p,r_p)\) forat least one surface point \(p\), where
\begin{equation}r_p = \max \{ r'>0 : B(p+r'n_p,r')\subset O\}\;.\label{eq:nch-continuous-radius}\end{equation}Note that in the mapping \(S\rightarrow \hbox{MAT}(O)\) so defined, twoor more surface points may map onto the same medial ball. But this isnot a problem, and in fact it is an unusual event; what is importantis that every medial ball is accounted for, so that the surface \(S\)can be reconstructed as the boundary of the union of balls\begin{equation}O = \cup \{B(p+r_p n_p,r_p) : p\in S \}\;.\label{eq:S-reconstruction-as union-of-balls}\end{equation}
The surface \(S\) can is approximated as the boundary of the Non-ConvexHull \(\hbox{NCH}({\cal P})\) defined as a half space of the NCH SignedDistance \(f(x)\). The Local Feature Size \(\hbox{LFS}(p)\) at a surfacepoint \(p\in S\) is usually defined as the distance from \(p\) to thenearest point in the Symmetric Medial Axis [Amenta et al. 2001; Dey 2007].A finite set \({\cal P}\subset S\) is defined as an \(\epsilon\)-sample of\(S\) if the distance from any point \(p\in S\) to its closest sample in\({\cal P}\) is at most \(\epsilon\,\hbox{LFS}(p)\). Several authors haveshown that for sufficiently small \(\epsilon\) the surface reconstructedas the boundary of \(\hbox{MAT}({\cal P})\) is a geometrically accurateapproximation of \(S\) [Amenta et al. 2001; Dey 2007], and ourexperiments validate these theoretical results.
Octrees and Isosurfaces
Since typically the NCH Signed Distance function has constant sign inlarge regions, one way to potentially reduce the computational cost ofthe naïve algorithm is to detect those areas and not to evaluatethe function there. Following a standard recursive space partitionapproach, we build an octree as a function of the point locations andthe orientation vectors, we evaluate the NCH Signed Distance functionat the vertices of the volumetric mesh, and use the Dual MarchingCubes (MC) algorithm [Schaefer and Warren 2005] to generate an outputpolygon mesh. Since the dual mesh of an octree is a conformingvolumetric polyhedral mesh, the polygon meshes produced by DMC areadaptive, but have no cracks. The results shown in figures 4 and 5have been computed using our implementation of DMC. Although theresults presented are very good, we regard them as preliminarywork. Due to lack of space, the details of this process as well asextensive experimental results will be presented in a future extendedpublication.
Conclusion
We have introduced an extremely simple algorithm to reconstructwatertight surfaces from finite sets of oriented points. Theformulation generalizes the Convex Hull in such a way that concavitiescan be represented and approximated. We have also proposedpreliminary strategies to reduce the computational cost by generatingadaptive polygon meshes and by subsampling. Despite its simplicity,the proposed algorithm produces high quality polygon meshescompetitive with those produced by state-of-the-art algorithms. Sincethe algorithm is massively paralellizable, and we plan to produce aGPU implementation in the near future.
Acknowledgements
This material is based upon work supported by the National ScienceFoundation under grants CCF-0729126, IIS-0808718, CCF-0915661, andIIP-1215308.
References
ALEXA, M., BEHR, J., COHEN-OR, D., FLEISHMAN, S., LEVIN, D., ANDSILVA, C. 2003. Computing and rendering point set surfaces. IEEETransactions on Visualization and Computer Graphics, 3–15.
ALLIEZ, P., COHEN-STEINER, D., TONG, Y., AND DESBRUN,M. 2007. Voronoi-based variational reconstruction of unoriented pointsets. In Proceedings of the fifth Eurographics symposium on Geometryprocessing, Eurographics Association, 39–48.
AMENTA, N., CHOI, S., AND KOLLURI, R. 2001. The Power Crust, Unions ofBalls, and the Medial Axis Transform. Computational Geometry Theoryand Applications 19, 2-3 (jul), 127– 153.
BERNARDINI, F., MITTLEMAN, J., RUSHMEIER, H., SILVA, C., AND TAUBIN,G. 1999. The Ball-Pivoting Algorithm for SurfaceReconstruction. IEEE Transactions on Visualization and ComputerGraphics 5, 4, 349–359.
BOISSONNAT, J., AND CAZALS, F. 2002. Smooth surface reconstruction vianatural neighbour interpolation of distance functions. ComputationalGeometry 22, 1, 185–203.
CALAKLI, F., AND TAUBIN, G. 2011. SSD: Smooth Signed DistanceSurface Reconstruction. Computer Graphics Forum 30,7. http://mesh.brown.edu/ssd.
DEY, T. 2007. Curve and surface reconstruction: algorithms withmathematical analysis. Cambridge University Press.
FLEISHMAN, S., COHEN-OR, D., AND SILVA, C. T. 2005. Robust movingleast-squares fitting with sharp features. ACM Transactions onGraphics 24 (July), 544–552.
HOPPE, H., DEROSE, T., DUCHAMP, T., MCDONALD, J., AND STUETZLE,W. 1992. Surface Reconstruction from Unorganized Points. InProceedings of ACM Siggraph, Citeseer.
KAZHDAN, M., BOLITHO, M., AND HOPPE, H. 2006. Poisson surfacereconstruction. In Proceedings of the fourth Eurographics symposiumon Geometry processing, Eurographics Association, 61–70.
LORENSEN, W., AND CLINE, H. 1987. Marching Cubes: A High Resolution 3dSurface Construction Algorithm. 163–169.
MANSON, J., PETROVA, G., AND SCHAEFER, S. 2008. Streaming surfacereconstruction using wavelets. In Computer Graphics Forum, vol. 27,Wiley Online Library, 1411–1420.
OHTAKE, Y., BELYAEV, A., ALEXA, M., TURK, G., AND SEIDEL,H. 2005. Multi-level partition of unity implicits. In ACM SIGGRAPH2005 Courses, ACM, 173.
SCHAEFER, S., AND WARREN, J. 2005. Dual Marching Cubes: primalcontouring of dual grids. In Computer Graphics Forum, vol. 24, WileyOnline Library, 195–201.