Feature Article
Iterative Consolidation of Unorganized Point Clouds Shengjun Liu, Kwan-Chung Chan, and Charlie C.L. Wang ■ Chinese University of Hong Kong
3
D scanners are now common; shapes obtained from them have become a major source of complex digital 3D models. Although optical scanners are the most economical and efficient way to obtain 3D digital models from real objects, they produce incomplete and noisy point clouds owing to occlusions and the scanners’ physical limitations. In the regions invisible to the cameras (deep cavities and bifurcations), the sample points don’t cover the scanned model’s surface. Most reconstruction methods lead to an imperfect surface Unorganized point clouds shape for the undersampled or obtained from 3D shape completely missed regions on the acquisition devices usually scanned point cloud of a realpresent noise, outliers, and world geometry (see Figure 1a). nonuniformities. A new For more information, see the framework consolidates “Related Work in Surface-Point unorganized points by Reconstruction” sidebar. interlaced downsampling and We’ve developed a framework for processing point clouds that upsampling. improves their quality and thus the reconstructed surfaces’ quality (see Figure 1b). The input to our framework is an unorganized point cloud that might contain outliers, noise, and nonuniformities. Based on the point positions, our framework evenly distributes points by inserting samples into sparse regions, using interlaced downsampling and upsampling and various operators. In particular, we use particle repulsion to ensure that the inserted sample points extrapolate the point set’s missing regions.
Overview of the Framework
Given an unorganized set of points P = {pj} ⊂ ℜ3 including outliers, noise, and nonuniformities, surface reconstruction from this data might cause 70
May/June 2012
significant misinterpretation, leading to an erroneous surface. Our framework aims to recover the structural information of P without losing the geometric details. It does this by inserting points in sparse regions to make the points evenly distributed, adding points to large missing regions, and removing outliers that are far from the upsampled surface. Such “massaging” of point clouds is called consolidation.5 Our framework comprises five steps—downsampling, outlier removal, repulsion, upsampling, and selection—which we iteratively apply to P. Our framework is inspired by an image-completion approach based on multiresolution techniques.6 That approach is based on the observation that a lower-resolution representation of an image contains stronger structural information and a higher-resolution representation contains more details. So, structural information is recovered at lower resolution. First, we downsample the points in P into k particles, which we redistribute on the surface defined by the samples in P. We perform this redistribution by iteratively applying a Weighted Locally Optimal Projection (WLOP) operator.7 A selection operator removes the outliers. We can estimate the oriented normals of particles by the cleaned particles.4,7 Our repulsion operator, which is based on the Algebraic Point Set Surface (APSS),8 extrapolates the surface by pushing particles into the missed regions. After that, we refine the redistributed particles into a smooth − point set surface using √3-like interpolatory refinement,9 which is an upsampling step. We selectively merge the newly generated sample points into P, treating the points in P that are far from the upsampled points as outliers and removing
Published by the IEEE Computer Society
0272-1716/12/$31.00 © 2012 IEEE
them. This produces a smooth surface that interpolates the redistributed particles. The downsampling, repulsion, upsampling, and selection steps are repeatedly applied to the point set. The iterations of this process end when either (P i–1−P i)/P i–1 < 20 percent or 10 iterations have occurred, with … being the number of points. Next, we obtain the orientation of the points in P by their closest particle in particle set X. We then obtain the consistently oriented normal vectors by orientation-aware principal component analysis. By itself, the repeated application of downsampling and upsampling to the whole set of point samples would waste much time on the regions that were processed in previous iterations. To reduce redundant computation, we use an adaptive approach. As Figure 2 shows, we downsample only the newly added points into particles while retaining the particles used in the previous iterations, and only upsample the newly added particles into new points. This approach speeds up the computation three to five times. Figure 3 illustrates the point-processing algorithm.
Original points (a)
POS
RBF
INT
OUR + POS
OUR + RBF
OUR + INT
Point-Processing Operators Here we discuss the technical details of the operators for downsampling and relaxation, mean-shiftbased outlier removal, repulsion, and upsampling.
Downsampling and Relaxation The WLOP operator performs downsampling and relaxation. For P i, we randomly select mi points to form set Xi. We call those points particles, and we select mi as mi = (2mP i)/P0, where m is a user parameter—we usually choose m = 1/20P 0 ∼ 1/5P 0. We then move these m particles to a new position, which is defined by two terms. The first attracts the particles to the given point set; the second repulses them away from other particles. For more details, see the “Position Update in Weighted Locally Optimal Projection” sidebar. We adopt a uniform support size h = 2Lavg for computing particle movement, where Lavg is the average distance between particles and their k nearest neighbors. So, the computation adapts to a model’s scale. We chose k = 20 to balance speed and robustness. When i ≠ 0, we move only the particles in Xi. All other particles in X are involved in the computation but have fixed positions. The choice of m is crucial. Using too few particles could easily damage the point set’s small features, whereas too many will slow the computation, relaxation, and repulsion of particles.7 So, we suggest m = 1/20P 0, for dense point clouds and m = 1/10P 0 for relatively sparse points.
OUR (b)
Figure 1. Comparing reconstruction methods. (a) The point cloud of an Inukshuk model (obtained from a 3D scanner) with incomplete sample points, along with surfaces reconstructed by a Poisson method (labeled POS),1 a radial-basis-function-based method (RBF),2 and an integrating-meshing method (INT).3 These methods are poor at regions with imperfect input samples. (b) The point cloud processed by our framework (labeled OUR), along with surfaces reconstructed from that cloud by the three methods. For the methods that need consistently oriented normal vectors, we generated them using Shengjun Liu and Charlie Wang’s method.4
Outlier Removal Although the WLOP operator can efficiently filter out noise by robust downsampling and relaxation, it can’t remove outliers that are far from a model’s real surface (see Figures 4a and 4b). These unwanted particles could affect the quality of the consolidation in the downstream operations and should be removed. Because these outliers are usually far from most of their neighbors, we employ a IEEE Computer Graphics and Applications
71
Feature Article
Figure 2. Our point-processing framework downsamples unorganized points (gray dots in the first column) into particles (yellow dots), redistributes them, and upsamples them into a dense point set (green dots). From the points generated by the upsampling, the framework selects and retains those in regions that lack samples in the original point set (red dots). In the next iteration, the framework downsamples the points retained from upsampling into new particles (the yellow dots in the second and third columns), redistributes them, and upsamples them into new points (green dots). The iterations repeat until only a few upsampled points are added.
mean-shift-based operator that detects and removes them during the consolidation’s first iteration. We iteratively shift each particle x in X to the average position of its k nearest neighbors (k = 20). Then, we compare the distance from the original position of x to the shifted mean position with the average distance between its shifted neighbors. If the difference is significant, we consider x an outlier. Figure 5 illustrates the mean-shift-based outlier removal algorithm. The average distance between shifted particles is computed locally, so it adapts to the nonuniform distribution of particles. In our experiments, we found that b = 3 and s = 3 give a good balance between speed and quality where b indicates the iterations number for the meanshift step and s is a user-defined parameter for the outlier detection. Figure 4c shows the results of removing outlier particles. Figure 4d shows that the mean-shift iteration converges quickly. So, 72
May/June 2012
using a larger b might not change the results of outlier identification or removal. Simply deleting the outlier particles that belong to the point set downsampled from P doesn’t remove the outliers from P because the outliers are only a small subset of the points in P. To remove them, we use a simple selection operation. After we obtain a set of cleaned particles and upsample them back into points in γ0 as the samples of a smooth surface (see Figure 4e), the points in P that are far from the surface represented by γ0 are considered outliers. So, for all points pi ∈ P, we remove a point p i when kp i – q jk > h (for all q j ∈ γ0), with q indicating the points in γ0. Figure 4f shows this operation’s results.
Particle Repulsion The WLOP operator evenly redistributes the particles along the surface defined by P. However,
Related Work in Surface-Point Reconstruction
R
esearchers have investigated reconstructing a surface from points for many years, and the reconstruction methods have become a standard for geometry creation. However, the points acquired by scanners are typically incomplete and highly nonuniform. Here we review methods for reconstructing surfaces from samples with inhomogeneous density or missing data. Yutaka Ohtake and his colleagues1 and Jonathan Carr and his colleagues2 exploit the extrapolation properties of radial basis functions to fill regions of sparse sampling. Vladimir Savchenko and Nikita Kojekine use control points to warp a given surface model toward that surface’s missing region.3 They then employ a fairing step along the hole’s boundary. This method requires manual interventions, and a prior model must be given in advance. Joan Verdera and his colleagues also use an implicit function to represent the surface.4 They model a partial differential equation for the smooth interpolation of a given hole on the basis of the normal vector field around it. Ulrich Clarenz and his colleagues use optimization to repair a surface.5 This process minimizes the integral of the squared mean curvature to yield a smooth surface. Tim Weyrich and his colleagues6 extend James Davis and his colleagues’ volumetric diffusion method7 to pointsampled models by replacing the distance estimation with a moving least-squares projection step. The previous methods can repair small deficiencies in the data but have difficulties when holes are complex or the object is missing large parts. Ravi Kolluri and his colleagues filter a Voronoi diagram to obtain a correct pole graph.8 To compute a watertight surface, they use global normalized cuts that smoothly complete large missing parts. The surface synthesis completes missing parts by integrating patches from a wellannotated shape database or a given example set.9,10 If no appropriate examples exist, the result might be poor and the process might fail. Alexander Hornung and Leif Kobbelt’s method requires defining a watertight voxel crust in which the unknown surface is supposed to lie.11 To complete the crust, they use flood-fill and dilation operators. Andrei Sharf and his colleagues evolve an explicit mesh in a scalar field guided by the local feature size in a coarse to fine manner to avoid local minima and capture details.12 Their method also requires a volumetric grid to evaluate the distance transformation, and the topological change must be tracked. The computational implementation can be quite intricate (especially the topology variation on two-manifold mesh surfaces). In other research, Sharf and his colleagues interactively reconstruct the surface using only the positions of raw scanned data. The user defines the general in and out orientations and interprets the data of automatically detected topologically unstable regions.13
Unlike the mesh surface reconstruction methods, our research (see the main article) focuses on improving the quality of an unorganized point cloud by adding sample points to change the scattered points’ uniformity and removing outliers. After iteratively consolidating the point cloud, we reconstruct a high-quality surface from those points by various methods.
References 1. Y. Ohtake, A. Belyaev, and H.-P. Seidel, “3D Scattered Data Interpolation and Approximation with Multilevel Compactly Supported RBFs,” Graphical Models, vol. 67, no. 3, 2005, pp. 150–165. 2. J.C. Carr et al., “Reconstruction and Representation of 3D Objects with Radial Basis Functions,” Proc. Siggraph, ACM Press, 2001, pp. 67–76. 3. V. Savchenko and N. Kojekine, “An Approach to Blend Surfaces,” Proc. Computer Graphics Int’l, Springer, 2002, pp. 139–150. 4. J. Verdera et al., “Inpainting Surface Holes,” Proc. Int’l Conf. Image Processing, IEEE CS Press, 2003, pp. 903–906. 5. U. Clarenz et al., “A Finite Element Method for Surface Restoration with Smooth Boundary Conditions,” Computer Aided Geometric Design, vol. 21, no. 5, 2004, pp. 427–455. 6. T. Weyrich et al., “Post-processing of Scanned 3D Surface Data,” Proc. Symp. Point-Based Graphics, Eurographics Assoc., 2004, pp. 85–94. 7. J. Davis et al., “Filling Holes in Complex Surfaces Using Volumetric Diffusion,” Proc. 1st Int’l Symp. 3D Data Processing, Visualization, and Transmission, IEEE CS Press, 2002, pp. 428–438. 8. R. Kolluri, J.R. Shewchuk, and J.F. O’Brien, “Spectral Surface Reconstruction from Noisy Point Clouds,” Proc. 2004 Eurographics/ACM Symp. Geometry Processing (SGP 04), ACM Press, 2004, pp. 11–21. 9. A. Sharf, M. Alexa, and D. Cohen-Or, “Context-Based Surface Completion,” ACM Trans. Graph., vol. 23, no. 3, 2004, pp. 878–887. 10. M. Pauly et al., “Example-Based 3D Scan Completion,” Proc. 3rd Eurographics Symp. Geometry Processing (SGP 05), Eurographics Assoc., 2005, pp. 23–32. 11. A. Hornung and L. Kobbelt, “Robust Reconstruction of Watertight 3D Models from Non-uniformly Sampled Point Clouds without Normal Information,” Proc. 4th Eurographics Symp. Geometry Processing (SGP 06), Eurographics Assoc., 2006, pp. 41–50. 12. A. Sharf et al., “Competing Fronts for Coarse-to-Fine Surface Reconstruction,” Proc. Eurographics, Eurographics Assoc., 2006, pp. 389–398. 13. A. Sharf et al., “Interactive Topology-Aware Surface Reconstruction,” ACM Trans. Graph., vol. 26, no. 3, 2007, article 43. IEEE Computer Graphics and Applications
73
Feature Article
Position Update in Weighted Locally Optimal Projection
G
iven an unorganized set of points P = {p j} ⊂ ℜ3, Weighted Locally Optimal Projection moves every particle xi ∈ Xk to a new position. The position update consists of two terms. The first term attracts the particle to the given point set by the weighted local density vj: v j = 1+
∑
pl ∈(P \{p j })
q ( p j − pl ).
The second term repulses the particles away from other particles by the density w p = 1+
∑
xq ∈(X \{x p })
q ( x p − xq ) ,
with q as a subscript of a particle, xq (distinguished from particle xp). The updated position of xi is xi =
∑
p j ∈P
pj
) v j xi − p j ∑ p ∈P θ ( xi − p j ) v j xi − p j w q θ ( xi − xq ) xi − xq x −x
∑
( x ∈(X \{x }) q
W=
∑
xq ∈(X \{xi })
w q q ( xi − xq
)
xi − xq , … 2
2
is the L 2-norm, and we adopt q (r ) = e −16r h as in Hui Huang and her colleagues’ method.1 q(r) is a rapidly decreasing smooth weight function with the support radius h defining the influenced neighborhood’s size. m ∈ [0, 0.5) and h are two user-selected parameters to tune the operator’s performance. During our tests, the default values m = 0.45 and h = 2Lavg worked well on all models, with Lavg being the average distance between particles and their k-nearest-neighbor particles. (We choose k = 20.)
θ ( xi − p j
Reference
j
+µ
where
i
i
q
)
W
,
1. H. Huang et al., “Consolidation of Unorganized Point Clouds for Surface Reconstruction,” ACM Trans. Graph., vol. 28, no. 5, 2009, pp. 1–7.
the movement of particles driven by WLOP stops at the boundary of large missing regions (see the missed region on the thighs in Figure 6). Other researchers have reported the same observation.7 To overcome this drawback, our repulsion operator moves the particles into the missing regions. These
regions’ normals are always nearly perpendicular to the scanner’s viewing direction. Smooth extension of the known regions that have been sampled could be a good heuristic to fill the missing regions. In other words, we must extrapolate the surface defined by the existing particles along its tangential direction.
P 0 ⇐ P and i ⇐ 0 Initialize a particle set X by downsampling all points of P into m particles repeat X ⇐ X ∪ Xi Repeatedly move the particles in Xi by the Weighted Locally Optimal Projection operator* if i = 0 then Remove the outlier particles from X0 by a mean-shift-based selection operation end if Estimate the particles’ orientation using Liu and Wang’s method4 or Hui Huang and her colleagues’ method7 Apply the repulsion operator based on Algebraic Point Set Surfaces to all particles in Xi Refine the points in Xi into a set of upsampled points γi if i = 0 then Remove the outliers in P according to γ0 end if Select the points of γi into a subset P i+1 P ⇐ P ∪ P i+1 and i ⇐ i + 1 Downsample all points of Pi into Xi with mi particles (mi = 2mP i/P 0) until the terminal condition is reached Estimate the consistently oriented normals on the sample points in P by Liu and Wang’s method or Huang and her colleagues’ method return P * The computation also uses the particles in (XnXi), although it moves only the particles in Xi.
Figure 3. The point-processing algorithm. Points or particles processed in previous iterations aren’t processed further, thus eliminating much computational redundancy. 74
May/June 2012
(a)
(b)
(c)
(f)
(g)
(d)
(e)
Figure 4. Outlier removal on an armadillo model. (a) The original point cloud had 20 percent noisy points randomly distributed over the range of 1/50 of the bounding box’s diagonal length. (b) Downsampling and relaxation based on Weight Locally Optimal Projection (WLOP) produced these particles. (c) Our meanshift-based outlier removal operator successfully identified and removed the outliers. (d) Line segments illustrate the mean-shifts on particles after one, three, and five iterations. The identified outliers are red. (e) Upsampling from the clean particles produced this result, which has less geometric detail. (f) Removing outliers in sample points and merging the upsampled points produced these results. (g) This image shows the particles to be used in the second iteration. Live particles (those that can be moved in the repulsion) sampled from the newly inserted points are red; static particles are blue.
After relaxation and outlier removal, the particles in X have become more uniform and less noisy. So, we can estimate the consistently oriented normal vectors on the particles using Shengjun Liu and Charlie Wang’s method4 or Hui Huang and her
colleagues’ method.7 Equipped with the oriented normals, we can define an APSS by all particles xi ∈ X and their normal vectors n xi . APSS is a kind of moving least-squares (MLS) surface. Instead of plane fitting, APSS directly fits higher-order IEEE Computer Graphics and Applications
75
Feature Article for all xi ∈ X do Initialize x′i ⇐ xi for j = 1 to b do Search in X for the k-nearest neighbors of x′i and let them be the subset K
x′i ⇐
1 K
∑
∀ki ∈K
∑
k i − x′i ∀ki ∈K
K
; {the average distance}
Figure 5. The mean-shift-based outlier removal algorithm. This algorithm can identify and remove the outliers successfully.
algebraic spheres. It yields more stable results than planar MLS at regions with high curvature, which have thin, sharp features. For more details about APSS fitting, see the related sidebar. After defining the APSS, we can move the particles along it in a repulsive manner. The movement has a tangential component and a projective component. We determine the tangential component of a particle xi by rotating along an axis rxi that passes through the center c xi of the algebraic sphere corresponding to xi. We derive the tangential component from the second term of WLOP (see the related sidebar) as
x p ∈( X \{x i })
(xi − x p )
∑
w p x p ∈( X \{x i })
w p
− xi ,
with = q ( x i − x p ) x i − x p and m as a parameter that balances the two terms of WLOP. It is described further in the “Position Update in Weighted Locally Optimal Projection” sidebar. Generally speaking, l is not perpendicular to n xi at x i. We thus compute the corresponding vector that’s orthogonal to n xi by l′ = l − (l ⋅ n xi ) n xi . The tangential component’s rotation axis is rxi =
( x i − c x )× l ′ , ( x i − c x )× l ′ i i
and we obtain the rotation angle by 76
May/June 2012
.
R rxi (ϖ )( x i − c xi ) + c xi − x i , where R rxi (ϖ ) is the rotation matrix10 around rxi . To improve the particle movement’s stability, we restrict the rotation angle by
After applying the tangential component, we apply the projective component on xi three times consecutively to retain the moved particles on the APSS defined by the particles with normal vectors:
return X
∑
2π x i − c xi
l′ Lavg π ϖ = min , , . 2π x i − c xi 2π x i − c xi 4
if x′i − xi > sd xi′ then xi is considered an outlier particle end if end for Remove all outlier particles from X
l=m
l′
We then define the tangential component as
k i ;{mean-shift step}
end for end for for all xi ∈ X do Search in X for the k-nearest neighbors of x′i and let them be the subset K
d xi′ ⇐
ϖ=
c xi + r ( x i )
x i − c xi − xi . x i − c xi
Applying this repulsion only once can’t fill large missing regions (see Figure 7). This is because the APSS defined by the particles before repulsion isn’t fully defined in the missing region (see the tessellation results in Figure 7d). So, we must fill the regions iteratively.
Upsampling and Selection Upsampling generates more sample points on the surface that interpolates the particles in Xi and the normal vectors on the particles. The upsampled point set γk should regularize the scattered samples and converge on a smooth surface that interpolates the particles and their normals. A good upsampling method that satisfies these requirements is Gaël Guennebaud and Markus Gross’s interpolatory refinement,8 which we use to generate γk from Xi. To fit our adaptive downsampling-upsampling strategy, we use only the live particles (those that can be moved in the repulsion) and the points upsampled from them as the centers to generate refined points. We prevent refinement at centers that are too close to their neighbors (for example, less than one-third of h). The upsampling stops when no refinement on any center is allowed. The repulsion and upsampling operators work together to generate the sample points on a smooth surface that extrapolates the WLOP-generated particles. Such extrapolation helps improve the uniformity of sample points in highly sparse regions. We select the points in γk to merge into the point set P k to form a new point set P k+1. In the next iteration of point processing, only the particles
(b)
(a)
(c) Figure 6. Results for the (a) Inukshuk model (b) without and (c) with the repulsion operator. Simply downsampling and upsampling the points without the operator couldn’t fill samples into the point cloud’s missing regions. Our framework, on the other hand, progressively improved the point cloud’s quality. Upsampled points are red.
(a)
(b)
(c)
(d)
Figure 7. How the repulsion operator works. (a) The given points for a human body, with large regions missing. (b) The redistributed particles after applying WLOP. (c) The redistributed particles after applying our repulsion operator. (d) The mesh surface tessellated from an Algebraic Point Set Surface, in which no triangle is generated in the undefined regions. One step of repulsion can fill some of the missing regions (for example, the arms). However, more iterations are needed to fill large missing regions (for example, the side of the head).
downsampled from these newly inserted points will be live—that is, moved during WLOP-based relaxation. So, simply adding all the points of γk to P k+1 would lead to too many particles in later iterations, which would significantly slow down the computation. Adding too few upsampled points would make too few live particles in the next iteration, which would also slow down the extrapolation generated by particle repulsion. Thus, we merge γk and P k with the criterion that for all q j ∈ γk, we exclude q j from P k+1 if there pnt exists pi ∈ Pk with q j − pi < 1 2 Lpnt avg , where Lavg is the average distance between all samples in P k and their k nearest neighbors. According to our experiments, using different values for 1 2 Lpnt avg doesn’t significantly affect the point set’s quality; it mainly affects the computation speed.
Results and Discussion We tested our framework on point clouds that had problems related to nonuniform sample points. Our experiments exhibited difficulties in reconstructing surfaces from the raw data directly in highly sparse areas, such as ■■ ■■
■■
■■
the legs of the Inukshuk in Figure 1, the parts between two fingers on a hand in Figure 8a, the armpit of the sculpture of a Japanese woman in Figure 8b, and the misscanned regions between the front and the back scans of the human body in Figure 9.
For illustration purposes, we rendered the unoriented raw scans using the normals computed by IEEE Computer Graphics and Applications
77
Feature Article Figure 8. Pointprocessing results. (a) A hand model. Previous methods have problems with the missed regions of the input point cloud between the fingers, but our framework can recover them. (b) A model of a sculpture of a Japanese woman. Our framework doesn’t change the region that was wellpresented in the input point cloud.
Original points + POS
(a)
(b)
OUR + POS
Original points + POS
Liu and Wang’s method.4 All the examples in this article take fewer than 10 iterations. Table 1 shows the computational statistics. For some point data that embeds small holes only, previous reconstruction methods can generate a surface without being significantly affected by the holes. In contrast, we aim to construct a high-quality point 78
May/June 2012
OUR + POS
set that can recover highly sparse regions while preserving geometric details. Toward that end, we use various surface reconstruction methods to reconstruct high-quality surfaces from the point clouds processed by our framework (see Figures 1 and 10). We also evaluated our framework’s robustness to noise. As Figure 11 shows, we applied it to the
Algebraic Point Set Surface Fitting
U
sing u = [u 0, …, u4] as the vector of scalar coefficients describing a general algebraic sphere in ℜ3, we evaluate the solution of algebraic-sphere fitting at a given point x ∈ ℜ3 with its normal nx by
u = arg min
u ,u≠0
1 W2
2
( x)Du .
W and D, respectively, are the 4n × 4n diagonal weight matrix and the 4n × 5 design matrix for the algebraic sphere fitting with normals; b is a 4n vector comprising of function values and normals of the fitting algebraic spheres at n points; and b is set by a large value for making the method insensitive to scale. This makes the method more stable, less prone to oscillations, and less sensitive to outliers. We solve for u by solving this linear-equation system:1 1
1
W 2 ( x) Du = W 2 ( x) b . The coefficient matrices have 4n rows, where n is the number of particles: w x ( ) i bw i ( x) , W ( x) = bw i ( x) w x b i( ) D=
1 xTi 0
eT0
0 eT1 0
eT2
2eT0 xi , 2eT1 xi 2eT2 xi
0 eT nx 0 i b= T e1 nxi eT nxi 2
.
The weight function w i ( x) = φ ( x − xi h) describes the weight of particle xi for the local evaluation of the Algebraic Point Set Surface at position x, with (1− r 2 )4 φ (r ) = 0
r
Iterative Consolidation of Unorganized Point Clouds Shengjun Liu, Kwan-Chung Chan, and Charlie C.L. Wang ■ Chinese University of Hong Kong
3
D scanners are now common; shapes obtained from them have become a major source of complex digital 3D models. Although optical scanners are the most economical and efficient way to obtain 3D digital models from real objects, they produce incomplete and noisy point clouds owing to occlusions and the scanners’ physical limitations. In the regions invisible to the cameras (deep cavities and bifurcations), the sample points don’t cover the scanned model’s surface. Most reconstruction methods lead to an imperfect surface Unorganized point clouds shape for the undersampled or obtained from 3D shape completely missed regions on the acquisition devices usually scanned point cloud of a realpresent noise, outliers, and world geometry (see Figure 1a). nonuniformities. A new For more information, see the framework consolidates “Related Work in Surface-Point unorganized points by Reconstruction” sidebar. interlaced downsampling and We’ve developed a framework for processing point clouds that upsampling. improves their quality and thus the reconstructed surfaces’ quality (see Figure 1b). The input to our framework is an unorganized point cloud that might contain outliers, noise, and nonuniformities. Based on the point positions, our framework evenly distributes points by inserting samples into sparse regions, using interlaced downsampling and upsampling and various operators. In particular, we use particle repulsion to ensure that the inserted sample points extrapolate the point set’s missing regions.
Overview of the Framework
Given an unorganized set of points P = {pj} ⊂ ℜ3 including outliers, noise, and nonuniformities, surface reconstruction from this data might cause 70
May/June 2012
significant misinterpretation, leading to an erroneous surface. Our framework aims to recover the structural information of P without losing the geometric details. It does this by inserting points in sparse regions to make the points evenly distributed, adding points to large missing regions, and removing outliers that are far from the upsampled surface. Such “massaging” of point clouds is called consolidation.5 Our framework comprises five steps—downsampling, outlier removal, repulsion, upsampling, and selection—which we iteratively apply to P. Our framework is inspired by an image-completion approach based on multiresolution techniques.6 That approach is based on the observation that a lower-resolution representation of an image contains stronger structural information and a higher-resolution representation contains more details. So, structural information is recovered at lower resolution. First, we downsample the points in P into k particles, which we redistribute on the surface defined by the samples in P. We perform this redistribution by iteratively applying a Weighted Locally Optimal Projection (WLOP) operator.7 A selection operator removes the outliers. We can estimate the oriented normals of particles by the cleaned particles.4,7 Our repulsion operator, which is based on the Algebraic Point Set Surface (APSS),8 extrapolates the surface by pushing particles into the missed regions. After that, we refine the redistributed particles into a smooth − point set surface using √3-like interpolatory refinement,9 which is an upsampling step. We selectively merge the newly generated sample points into P, treating the points in P that are far from the upsampled points as outliers and removing
Published by the IEEE Computer Society
0272-1716/12/$31.00 © 2012 IEEE
them. This produces a smooth surface that interpolates the redistributed particles. The downsampling, repulsion, upsampling, and selection steps are repeatedly applied to the point set. The iterations of this process end when either (P i–1−P i)/P i–1 < 20 percent or 10 iterations have occurred, with … being the number of points. Next, we obtain the orientation of the points in P by their closest particle in particle set X. We then obtain the consistently oriented normal vectors by orientation-aware principal component analysis. By itself, the repeated application of downsampling and upsampling to the whole set of point samples would waste much time on the regions that were processed in previous iterations. To reduce redundant computation, we use an adaptive approach. As Figure 2 shows, we downsample only the newly added points into particles while retaining the particles used in the previous iterations, and only upsample the newly added particles into new points. This approach speeds up the computation three to five times. Figure 3 illustrates the point-processing algorithm.
Original points (a)
POS
RBF
INT
OUR + POS
OUR + RBF
OUR + INT
Point-Processing Operators Here we discuss the technical details of the operators for downsampling and relaxation, mean-shiftbased outlier removal, repulsion, and upsampling.
Downsampling and Relaxation The WLOP operator performs downsampling and relaxation. For P i, we randomly select mi points to form set Xi. We call those points particles, and we select mi as mi = (2mP i)/P0, where m is a user parameter—we usually choose m = 1/20P 0 ∼ 1/5P 0. We then move these m particles to a new position, which is defined by two terms. The first attracts the particles to the given point set; the second repulses them away from other particles. For more details, see the “Position Update in Weighted Locally Optimal Projection” sidebar. We adopt a uniform support size h = 2Lavg for computing particle movement, where Lavg is the average distance between particles and their k nearest neighbors. So, the computation adapts to a model’s scale. We chose k = 20 to balance speed and robustness. When i ≠ 0, we move only the particles in Xi. All other particles in X are involved in the computation but have fixed positions. The choice of m is crucial. Using too few particles could easily damage the point set’s small features, whereas too many will slow the computation, relaxation, and repulsion of particles.7 So, we suggest m = 1/20P 0, for dense point clouds and m = 1/10P 0 for relatively sparse points.
OUR (b)
Figure 1. Comparing reconstruction methods. (a) The point cloud of an Inukshuk model (obtained from a 3D scanner) with incomplete sample points, along with surfaces reconstructed by a Poisson method (labeled POS),1 a radial-basis-function-based method (RBF),2 and an integrating-meshing method (INT).3 These methods are poor at regions with imperfect input samples. (b) The point cloud processed by our framework (labeled OUR), along with surfaces reconstructed from that cloud by the three methods. For the methods that need consistently oriented normal vectors, we generated them using Shengjun Liu and Charlie Wang’s method.4
Outlier Removal Although the WLOP operator can efficiently filter out noise by robust downsampling and relaxation, it can’t remove outliers that are far from a model’s real surface (see Figures 4a and 4b). These unwanted particles could affect the quality of the consolidation in the downstream operations and should be removed. Because these outliers are usually far from most of their neighbors, we employ a IEEE Computer Graphics and Applications
71
Feature Article
Figure 2. Our point-processing framework downsamples unorganized points (gray dots in the first column) into particles (yellow dots), redistributes them, and upsamples them into a dense point set (green dots). From the points generated by the upsampling, the framework selects and retains those in regions that lack samples in the original point set (red dots). In the next iteration, the framework downsamples the points retained from upsampling into new particles (the yellow dots in the second and third columns), redistributes them, and upsamples them into new points (green dots). The iterations repeat until only a few upsampled points are added.
mean-shift-based operator that detects and removes them during the consolidation’s first iteration. We iteratively shift each particle x in X to the average position of its k nearest neighbors (k = 20). Then, we compare the distance from the original position of x to the shifted mean position with the average distance between its shifted neighbors. If the difference is significant, we consider x an outlier. Figure 5 illustrates the mean-shift-based outlier removal algorithm. The average distance between shifted particles is computed locally, so it adapts to the nonuniform distribution of particles. In our experiments, we found that b = 3 and s = 3 give a good balance between speed and quality where b indicates the iterations number for the meanshift step and s is a user-defined parameter for the outlier detection. Figure 4c shows the results of removing outlier particles. Figure 4d shows that the mean-shift iteration converges quickly. So, 72
May/June 2012
using a larger b might not change the results of outlier identification or removal. Simply deleting the outlier particles that belong to the point set downsampled from P doesn’t remove the outliers from P because the outliers are only a small subset of the points in P. To remove them, we use a simple selection operation. After we obtain a set of cleaned particles and upsample them back into points in γ0 as the samples of a smooth surface (see Figure 4e), the points in P that are far from the surface represented by γ0 are considered outliers. So, for all points pi ∈ P, we remove a point p i when kp i – q jk > h (for all q j ∈ γ0), with q indicating the points in γ0. Figure 4f shows this operation’s results.
Particle Repulsion The WLOP operator evenly redistributes the particles along the surface defined by P. However,
Related Work in Surface-Point Reconstruction
R
esearchers have investigated reconstructing a surface from points for many years, and the reconstruction methods have become a standard for geometry creation. However, the points acquired by scanners are typically incomplete and highly nonuniform. Here we review methods for reconstructing surfaces from samples with inhomogeneous density or missing data. Yutaka Ohtake and his colleagues1 and Jonathan Carr and his colleagues2 exploit the extrapolation properties of radial basis functions to fill regions of sparse sampling. Vladimir Savchenko and Nikita Kojekine use control points to warp a given surface model toward that surface’s missing region.3 They then employ a fairing step along the hole’s boundary. This method requires manual interventions, and a prior model must be given in advance. Joan Verdera and his colleagues also use an implicit function to represent the surface.4 They model a partial differential equation for the smooth interpolation of a given hole on the basis of the normal vector field around it. Ulrich Clarenz and his colleagues use optimization to repair a surface.5 This process minimizes the integral of the squared mean curvature to yield a smooth surface. Tim Weyrich and his colleagues6 extend James Davis and his colleagues’ volumetric diffusion method7 to pointsampled models by replacing the distance estimation with a moving least-squares projection step. The previous methods can repair small deficiencies in the data but have difficulties when holes are complex or the object is missing large parts. Ravi Kolluri and his colleagues filter a Voronoi diagram to obtain a correct pole graph.8 To compute a watertight surface, they use global normalized cuts that smoothly complete large missing parts. The surface synthesis completes missing parts by integrating patches from a wellannotated shape database or a given example set.9,10 If no appropriate examples exist, the result might be poor and the process might fail. Alexander Hornung and Leif Kobbelt’s method requires defining a watertight voxel crust in which the unknown surface is supposed to lie.11 To complete the crust, they use flood-fill and dilation operators. Andrei Sharf and his colleagues evolve an explicit mesh in a scalar field guided by the local feature size in a coarse to fine manner to avoid local minima and capture details.12 Their method also requires a volumetric grid to evaluate the distance transformation, and the topological change must be tracked. The computational implementation can be quite intricate (especially the topology variation on two-manifold mesh surfaces). In other research, Sharf and his colleagues interactively reconstruct the surface using only the positions of raw scanned data. The user defines the general in and out orientations and interprets the data of automatically detected topologically unstable regions.13
Unlike the mesh surface reconstruction methods, our research (see the main article) focuses on improving the quality of an unorganized point cloud by adding sample points to change the scattered points’ uniformity and removing outliers. After iteratively consolidating the point cloud, we reconstruct a high-quality surface from those points by various methods.
References 1. Y. Ohtake, A. Belyaev, and H.-P. Seidel, “3D Scattered Data Interpolation and Approximation with Multilevel Compactly Supported RBFs,” Graphical Models, vol. 67, no. 3, 2005, pp. 150–165. 2. J.C. Carr et al., “Reconstruction and Representation of 3D Objects with Radial Basis Functions,” Proc. Siggraph, ACM Press, 2001, pp. 67–76. 3. V. Savchenko and N. Kojekine, “An Approach to Blend Surfaces,” Proc. Computer Graphics Int’l, Springer, 2002, pp. 139–150. 4. J. Verdera et al., “Inpainting Surface Holes,” Proc. Int’l Conf. Image Processing, IEEE CS Press, 2003, pp. 903–906. 5. U. Clarenz et al., “A Finite Element Method for Surface Restoration with Smooth Boundary Conditions,” Computer Aided Geometric Design, vol. 21, no. 5, 2004, pp. 427–455. 6. T. Weyrich et al., “Post-processing of Scanned 3D Surface Data,” Proc. Symp. Point-Based Graphics, Eurographics Assoc., 2004, pp. 85–94. 7. J. Davis et al., “Filling Holes in Complex Surfaces Using Volumetric Diffusion,” Proc. 1st Int’l Symp. 3D Data Processing, Visualization, and Transmission, IEEE CS Press, 2002, pp. 428–438. 8. R. Kolluri, J.R. Shewchuk, and J.F. O’Brien, “Spectral Surface Reconstruction from Noisy Point Clouds,” Proc. 2004 Eurographics/ACM Symp. Geometry Processing (SGP 04), ACM Press, 2004, pp. 11–21. 9. A. Sharf, M. Alexa, and D. Cohen-Or, “Context-Based Surface Completion,” ACM Trans. Graph., vol. 23, no. 3, 2004, pp. 878–887. 10. M. Pauly et al., “Example-Based 3D Scan Completion,” Proc. 3rd Eurographics Symp. Geometry Processing (SGP 05), Eurographics Assoc., 2005, pp. 23–32. 11. A. Hornung and L. Kobbelt, “Robust Reconstruction of Watertight 3D Models from Non-uniformly Sampled Point Clouds without Normal Information,” Proc. 4th Eurographics Symp. Geometry Processing (SGP 06), Eurographics Assoc., 2006, pp. 41–50. 12. A. Sharf et al., “Competing Fronts for Coarse-to-Fine Surface Reconstruction,” Proc. Eurographics, Eurographics Assoc., 2006, pp. 389–398. 13. A. Sharf et al., “Interactive Topology-Aware Surface Reconstruction,” ACM Trans. Graph., vol. 26, no. 3, 2007, article 43. IEEE Computer Graphics and Applications
73
Feature Article
Position Update in Weighted Locally Optimal Projection
G
iven an unorganized set of points P = {p j} ⊂ ℜ3, Weighted Locally Optimal Projection moves every particle xi ∈ Xk to a new position. The position update consists of two terms. The first term attracts the particle to the given point set by the weighted local density vj: v j = 1+
∑
pl ∈(P \{p j })
q ( p j − pl ).
The second term repulses the particles away from other particles by the density w p = 1+
∑
xq ∈(X \{x p })
q ( x p − xq ) ,
with q as a subscript of a particle, xq (distinguished from particle xp). The updated position of xi is xi =
∑
p j ∈P
pj
) v j xi − p j ∑ p ∈P θ ( xi − p j ) v j xi − p j w q θ ( xi − xq ) xi − xq x −x
∑
( x ∈(X \{x }) q
W=
∑
xq ∈(X \{xi })
w q q ( xi − xq
)
xi − xq , … 2
2
is the L 2-norm, and we adopt q (r ) = e −16r h as in Hui Huang and her colleagues’ method.1 q(r) is a rapidly decreasing smooth weight function with the support radius h defining the influenced neighborhood’s size. m ∈ [0, 0.5) and h are two user-selected parameters to tune the operator’s performance. During our tests, the default values m = 0.45 and h = 2Lavg worked well on all models, with Lavg being the average distance between particles and their k-nearest-neighbor particles. (We choose k = 20.)
θ ( xi − p j
Reference
j
+µ
where
i
i
q
)
W
,
1. H. Huang et al., “Consolidation of Unorganized Point Clouds for Surface Reconstruction,” ACM Trans. Graph., vol. 28, no. 5, 2009, pp. 1–7.
the movement of particles driven by WLOP stops at the boundary of large missing regions (see the missed region on the thighs in Figure 6). Other researchers have reported the same observation.7 To overcome this drawback, our repulsion operator moves the particles into the missing regions. These
regions’ normals are always nearly perpendicular to the scanner’s viewing direction. Smooth extension of the known regions that have been sampled could be a good heuristic to fill the missing regions. In other words, we must extrapolate the surface defined by the existing particles along its tangential direction.
P 0 ⇐ P and i ⇐ 0 Initialize a particle set X by downsampling all points of P into m particles repeat X ⇐ X ∪ Xi Repeatedly move the particles in Xi by the Weighted Locally Optimal Projection operator* if i = 0 then Remove the outlier particles from X0 by a mean-shift-based selection operation end if Estimate the particles’ orientation using Liu and Wang’s method4 or Hui Huang and her colleagues’ method7 Apply the repulsion operator based on Algebraic Point Set Surfaces to all particles in Xi Refine the points in Xi into a set of upsampled points γi if i = 0 then Remove the outliers in P according to γ0 end if Select the points of γi into a subset P i+1 P ⇐ P ∪ P i+1 and i ⇐ i + 1 Downsample all points of Pi into Xi with mi particles (mi = 2mP i/P 0) until the terminal condition is reached Estimate the consistently oriented normals on the sample points in P by Liu and Wang’s method or Huang and her colleagues’ method return P * The computation also uses the particles in (XnXi), although it moves only the particles in Xi.
Figure 3. The point-processing algorithm. Points or particles processed in previous iterations aren’t processed further, thus eliminating much computational redundancy. 74
May/June 2012
(a)
(b)
(c)
(f)
(g)
(d)
(e)
Figure 4. Outlier removal on an armadillo model. (a) The original point cloud had 20 percent noisy points randomly distributed over the range of 1/50 of the bounding box’s diagonal length. (b) Downsampling and relaxation based on Weight Locally Optimal Projection (WLOP) produced these particles. (c) Our meanshift-based outlier removal operator successfully identified and removed the outliers. (d) Line segments illustrate the mean-shifts on particles after one, three, and five iterations. The identified outliers are red. (e) Upsampling from the clean particles produced this result, which has less geometric detail. (f) Removing outliers in sample points and merging the upsampled points produced these results. (g) This image shows the particles to be used in the second iteration. Live particles (those that can be moved in the repulsion) sampled from the newly inserted points are red; static particles are blue.
After relaxation and outlier removal, the particles in X have become more uniform and less noisy. So, we can estimate the consistently oriented normal vectors on the particles using Shengjun Liu and Charlie Wang’s method4 or Hui Huang and her
colleagues’ method.7 Equipped with the oriented normals, we can define an APSS by all particles xi ∈ X and their normal vectors n xi . APSS is a kind of moving least-squares (MLS) surface. Instead of plane fitting, APSS directly fits higher-order IEEE Computer Graphics and Applications
75
Feature Article for all xi ∈ X do Initialize x′i ⇐ xi for j = 1 to b do Search in X for the k-nearest neighbors of x′i and let them be the subset K
x′i ⇐
1 K
∑
∀ki ∈K
∑
k i − x′i ∀ki ∈K
K
; {the average distance}
Figure 5. The mean-shift-based outlier removal algorithm. This algorithm can identify and remove the outliers successfully.
algebraic spheres. It yields more stable results than planar MLS at regions with high curvature, which have thin, sharp features. For more details about APSS fitting, see the related sidebar. After defining the APSS, we can move the particles along it in a repulsive manner. The movement has a tangential component and a projective component. We determine the tangential component of a particle xi by rotating along an axis rxi that passes through the center c xi of the algebraic sphere corresponding to xi. We derive the tangential component from the second term of WLOP (see the related sidebar) as
x p ∈( X \{x i })
(xi − x p )
∑
w p x p ∈( X \{x i })
w p
− xi ,
with = q ( x i − x p ) x i − x p and m as a parameter that balances the two terms of WLOP. It is described further in the “Position Update in Weighted Locally Optimal Projection” sidebar. Generally speaking, l is not perpendicular to n xi at x i. We thus compute the corresponding vector that’s orthogonal to n xi by l′ = l − (l ⋅ n xi ) n xi . The tangential component’s rotation axis is rxi =
( x i − c x )× l ′ , ( x i − c x )× l ′ i i
and we obtain the rotation angle by 76
May/June 2012
.
R rxi (ϖ )( x i − c xi ) + c xi − x i , where R rxi (ϖ ) is the rotation matrix10 around rxi . To improve the particle movement’s stability, we restrict the rotation angle by
After applying the tangential component, we apply the projective component on xi three times consecutively to retain the moved particles on the APSS defined by the particles with normal vectors:
return X
∑
2π x i − c xi
l′ Lavg π ϖ = min , , . 2π x i − c xi 2π x i − c xi 4
if x′i − xi > sd xi′ then xi is considered an outlier particle end if end for Remove all outlier particles from X
l=m
l′
We then define the tangential component as
k i ;{mean-shift step}
end for end for for all xi ∈ X do Search in X for the k-nearest neighbors of x′i and let them be the subset K
d xi′ ⇐
ϖ=
c xi + r ( x i )
x i − c xi − xi . x i − c xi
Applying this repulsion only once can’t fill large missing regions (see Figure 7). This is because the APSS defined by the particles before repulsion isn’t fully defined in the missing region (see the tessellation results in Figure 7d). So, we must fill the regions iteratively.
Upsampling and Selection Upsampling generates more sample points on the surface that interpolates the particles in Xi and the normal vectors on the particles. The upsampled point set γk should regularize the scattered samples and converge on a smooth surface that interpolates the particles and their normals. A good upsampling method that satisfies these requirements is Gaël Guennebaud and Markus Gross’s interpolatory refinement,8 which we use to generate γk from Xi. To fit our adaptive downsampling-upsampling strategy, we use only the live particles (those that can be moved in the repulsion) and the points upsampled from them as the centers to generate refined points. We prevent refinement at centers that are too close to their neighbors (for example, less than one-third of h). The upsampling stops when no refinement on any center is allowed. The repulsion and upsampling operators work together to generate the sample points on a smooth surface that extrapolates the WLOP-generated particles. Such extrapolation helps improve the uniformity of sample points in highly sparse regions. We select the points in γk to merge into the point set P k to form a new point set P k+1. In the next iteration of point processing, only the particles
(b)
(a)
(c) Figure 6. Results for the (a) Inukshuk model (b) without and (c) with the repulsion operator. Simply downsampling and upsampling the points without the operator couldn’t fill samples into the point cloud’s missing regions. Our framework, on the other hand, progressively improved the point cloud’s quality. Upsampled points are red.
(a)
(b)
(c)
(d)
Figure 7. How the repulsion operator works. (a) The given points for a human body, with large regions missing. (b) The redistributed particles after applying WLOP. (c) The redistributed particles after applying our repulsion operator. (d) The mesh surface tessellated from an Algebraic Point Set Surface, in which no triangle is generated in the undefined regions. One step of repulsion can fill some of the missing regions (for example, the arms). However, more iterations are needed to fill large missing regions (for example, the side of the head).
downsampled from these newly inserted points will be live—that is, moved during WLOP-based relaxation. So, simply adding all the points of γk to P k+1 would lead to too many particles in later iterations, which would significantly slow down the computation. Adding too few upsampled points would make too few live particles in the next iteration, which would also slow down the extrapolation generated by particle repulsion. Thus, we merge γk and P k with the criterion that for all q j ∈ γk, we exclude q j from P k+1 if there pnt exists pi ∈ Pk with q j − pi < 1 2 Lpnt avg , where Lavg is the average distance between all samples in P k and their k nearest neighbors. According to our experiments, using different values for 1 2 Lpnt avg doesn’t significantly affect the point set’s quality; it mainly affects the computation speed.
Results and Discussion We tested our framework on point clouds that had problems related to nonuniform sample points. Our experiments exhibited difficulties in reconstructing surfaces from the raw data directly in highly sparse areas, such as ■■ ■■
■■
■■
the legs of the Inukshuk in Figure 1, the parts between two fingers on a hand in Figure 8a, the armpit of the sculpture of a Japanese woman in Figure 8b, and the misscanned regions between the front and the back scans of the human body in Figure 9.
For illustration purposes, we rendered the unoriented raw scans using the normals computed by IEEE Computer Graphics and Applications
77
Feature Article Figure 8. Pointprocessing results. (a) A hand model. Previous methods have problems with the missed regions of the input point cloud between the fingers, but our framework can recover them. (b) A model of a sculpture of a Japanese woman. Our framework doesn’t change the region that was wellpresented in the input point cloud.
Original points + POS
(a)
(b)
OUR + POS
Original points + POS
Liu and Wang’s method.4 All the examples in this article take fewer than 10 iterations. Table 1 shows the computational statistics. For some point data that embeds small holes only, previous reconstruction methods can generate a surface without being significantly affected by the holes. In contrast, we aim to construct a high-quality point 78
May/June 2012
OUR + POS
set that can recover highly sparse regions while preserving geometric details. Toward that end, we use various surface reconstruction methods to reconstruct high-quality surfaces from the point clouds processed by our framework (see Figures 1 and 10). We also evaluated our framework’s robustness to noise. As Figure 11 shows, we applied it to the
Algebraic Point Set Surface Fitting
U
sing u = [u 0, …, u4] as the vector of scalar coefficients describing a general algebraic sphere in ℜ3, we evaluate the solution of algebraic-sphere fitting at a given point x ∈ ℜ3 with its normal nx by
u = arg min
u ,u≠0
1 W2
2
( x)Du .
W and D, respectively, are the 4n × 4n diagonal weight matrix and the 4n × 5 design matrix for the algebraic sphere fitting with normals; b is a 4n vector comprising of function values and normals of the fitting algebraic spheres at n points; and b is set by a large value for making the method insensitive to scale. This makes the method more stable, less prone to oscillations, and less sensitive to outliers. We solve for u by solving this linear-equation system:1 1
1
W 2 ( x) Du = W 2 ( x) b . The coefficient matrices have 4n rows, where n is the number of particles: w x ( ) i bw i ( x) , W ( x) = bw i ( x) w x b i( ) D=
1 xTi 0
eT0
0 eT1 0
eT2
2eT0 xi , 2eT1 xi 2eT2 xi
0 eT nx 0 i b= T e1 nxi eT nxi 2
.
The weight function w i ( x) = φ ( x − xi h) describes the weight of particle xi for the local evaluation of the Algebraic Point Set Surface at position x, with (1− r 2 )4 φ (r ) = 0
r
