LSMAT Least Squares Medial Axis Transform
Daniel Rebain, Baptiste Angles, Julien Valentin, Nicholas Vining, Jiju Peethambaran, Shahram Izadi, Andrea Tagliasacchi
LLSMATLeast Squares Medial Axis Transform
Daniel Rebain , , Baptiste Angles , , Julien Valentin , Nicholas Vining , Jiju Peethambaran , Shahram Izadi , Andrea Tagliasacchi , , Department of Computer Science, University of Victoria, Victoria, BC, Canada Google LLC, Mountain View, CA, USA David R. Cheriton School of Computer Science, University of Waterloo, Waterloo, ON, Canada
Figure 1: The LSMAT is a novel efficient least-squares formulation of the medial axis transform that operates on unorganized oriented pointsets. (a, b) Regardless of input noise, the resulting medial representation is stable . (c, d) Its least-squares nature allows it to operate in thepresence of heavy noise where most approaches would fail. We visualize the oriented point cloud with oriented splats, draw the union ofmedial spheres in light red, and their corresponding centers in dark red.
Abstract
The medial axis transform has applications in numerous fields including visualization, computer graphics, and computer vision.Unfortunately, traditional medial axis transformations are usually brittle in the presence of outliers, perturbations and/ornoise along the boundary of objects. To overcome this limitation, we introduce a new formulation of the medial axis transformwhich is naturally robust in the presence of these artifacts. Unlike previous work which has approached the medial axis from acomputational geometry angle, we consider it from a numerical optimization perspective. In this work, we follow the definition ofthe medial axis transform as “the set of maximally inscribed spheres". We show how this definition can be formulated as a leastsquares relaxation where the transform is obtained by minimizing a continuous optimization problem. The proposed approach isinherently parallelizable by performing independant optimization of each sphere using Gauss-Newton, and its least-squares formallows it to be significantly more robust compared to traditional computational geometry approaches. Extensive experiments on2D and 3D objects demonstrate that our method provides superior results to the state of the art on both synthetic and real-data.
Keywords: [medial axis transform] [optimization] [least squares]
CCS Concepts • Computing methodologies → Point-based models;
Volumetric models;
Shape analysis;
1. Introduction
A medial representation of an object encodes its solid geometry asthe union of a collection of spheres of different radii and origins;see Figure 2. While volumetric or surface mesh representations aremore commonly used in computer graphics and computer vision, since its introduction by Blum et al. [Blu67], medial representationshave found applications in many 2D/3D geometric problems such asanimation [BP07], fabrication [MHR ∗ The definitive version is available at wileyonlinelibrary.com a r X i v : . [ c s . G R ] O c t ebain et al. / LSMAT: Least Squares Medial Axis Transform maximally inscribed spheres normal flow shock graph bisectors / equidistance bi-tangency / local symmetry Figure 2: Visualization of four alternative definitions of medial axis; base image courtesy of [TDS ∗ medial skeleton . Stability.
A common pitfall of medial axis methods is their sensitiv-ity to noise: in a traditional medial representation, new branches ofthe medial skeleton may form at all negative local curvature extrema.If our shape is represented via a piecewise linear boundary, such as apolygonal mesh in 3D, spurious branches form when the surface ex-hibits even minor levels of noise; see [TDS ∗
16, Fig.9]. These issuesare typically resolved by resorting to postprocessing, in which setsof spheres are filtered out according to various engineered criteria.We present a new algorithm for computing the medial axis transformof an oriented point set which is naturally capable of handling noise,outliers, and other artifacts that create characteristic problems fortraditional methods.
Definitions.
As summarized in Figure 2 and [TDS ∗ normal flow variant leadsto the commonly employed voxel thinning algorithms [SBdB16],while the equidistance definition leads to techniques leveragingVoronoi diagrams [BA92]. In this paper, we build over what is likelythe most well known definition of the medial axis transform: Definition 1.1.
The Medial Axis Transform MAT ( O ) of O is the set of centers M and corresponding radii R ofall maximally inscribed circles/spheres in O . While previous work such as Ma et al. [MBC12] has proposed meth-ods to construct medial axis representations using this definition, itis interesting to note that they have treated the problem from a com-binatorial geometry standpoint. Instead, we consider the medial axistransform from a numerical geometry perspective, where the medialaxis is given by the solution of an optimization problem. We achievethis by expressing the concepts of maximality and inscription fromDef. 1.1 in a least squares form. The robustness of the approach toimperfections arises from the fact that least squares optimizationattempts to find an approximate, rather than exact, solution to thegiven problem. We are motivated in our approach by consideringa least-squares problem as a maximum likelihood estimate of a function in the presence of noise obeying a Gaussian probabilitydistribution.
Method outline.
Our method takes an oriented point cloud as input,and produces an unconnected medial point cloud as output by min-imizing a combination of a maximality energy, and an inscriptionenergy. Our key technical challenge is in formulating these energiescorrectly; our maximality energy is designed to have a constantmagnitude, ensuring that the optimization energy of each medialsphere is constant regardless of its radius, and our inscription energyis based around a locally supported approximation of the signeddistance function of the point cloud. In order to prevent spheres fromsliding towards local maxima of local shape thickness, we introducea pinning constraint as a quadratic barrier energy. The resultingoptimization problem is quadratic, with differentiable but non-linearenergy terms, and we solve it with an iterative Gauss-Newton solver.
Contributions and Evaluation.
Our main contributions are a novelinterpretation of a problem that is classically solved by standardcomputational geometry as a numerical geometry problem, anda novel algorithm for computing the medial axis transform of anoriented point set that inherently handles imperfections in the input.We present a number of qualitative results throughout the paper forboth 2D and 3D oriented point clouds. Finally, we also present side-by-side quantitative evaluations of the proposed algorithm againstseveral state-of-the-art methods.
2. Related Work
Many techniques exist to compute the medial axis transformation,as detailed in a recent survey [TDS ∗ surface instability of medial axis via postprocessing. Assuming the surface is sampled by a sufficiently dense set of points,the Voronoi diagram can be used to compute the medial axis trans-form with relative ease. For a closed boundary in 2D, any interiorVoronoi vertex, as well as Voronoi edges connecting interior vertices,approximate the medial axis with a convergence guarantee [BA92].
The definitive version is available at wileyonlinelibrary.com ebain et al. / LSMAT: Least Squares Medial Axis Transform unfiltered
Figure 3: Effectiveness of standard filtering Voronoi methods withhigh level of noise. As expected, neither angle nor distance methodswere capable of filtering out all noise for any selected thresholdvalue. Leveraging global information, the scale axis is able to dealwith noise, but with a significant computational burden. Further,note that we provide these technique the ground truth inside/outsidelabeling, which is not available in our context.Unfortunately, Amenta et al. [AB99] has shown how this prop-erty does not hold in 3D due to sliver tetrahedra . Thankfully, aswe increase the sampling rate, the Voronoi poles (see definitionin [AB99]) do converge to the medial axis, allowing the use of3D Voronoi diagrams for the task at hand. These methods suffertwo fundamental shortcomings: (cid:172) they are global and optimize theentire set of centers at the same time, and (cid:173) the extracted axis isheavily susceptible to even minor levels of noise; see Figure 3.
Sphere-at-a-point methods.
Leveraging an oriented sampling,Ma et al. [MBC12] proposed an alternative way to compute thetransform by marrying the maximally inscribed definition of Fig-ure 2(a) to the bi-tangency of Figure 2(d). Unlike Voronoi methods,which consider the entire point set at once, their method can com-pute a maximal sphere at a point in isolation , leading to extremelyefficient GPU implementations [JKT13]. These methods representthe most efficient way of computing the medial axis, but, similarlyto Voronoi methods, they suffer stability issues; see Figure 4. An-other method for computing a local approximation to the axis wasproposed by Shapira et al. [SSCO08] via casting rays in a coneoriented along the anti-normal. The distance between the point andthe intersection with the surface is aggregated by a robust function(e.g. median) to estimate the shape diameter function . While thisapproximation has found widespread use as a shape descriptor, theradius estimate suffer of bias , and the algorithm does not generalizeto point clouds.
Shape approximation methods.
Recently, a new class of tech-niques has been proposed which attempt to approximate water-tight surfaces via
Sphere Meshes – linearly swept spherical primi-tives [TGB13,TGBE16]. This is achieved via local mesh decimationrelying on iterative edge collapses, where spherical quadrics are em-ployed in place of the traditional quadric metrics [GH97]. In manycases the produced model resembles a medial axis. When executedon 3D data the result can contain tetrahedra, however, the medialaxis of a 3D shape is known to consist only of points, curves, andsurfaces. Addressing these concerns, the
Medial Meshes work bySun et al. [SCYW16] extended sphere meshes to decimate a medialaxis mesh, and discard unstable branches. Marrying sphere meshesto medial meshes, Li et al. [LWS ∗
15] followed up this work andproposed QMAT, a more computationally efficient version based onspherical quadrics. While these techniques, in juxtaposition to ourlocal method, are global, they can only cope with minor levels ofnoise: “with very noisy input, however, the simplified medial mesh isnot a stable representation” ; see [SCYW16, Fig.9]. The follow-upwork by Li et al. [LWS ∗
15] performs slightly better as it can opti-mize for sphere centers, but our algorithm can still cope with noisethat is one order of magnitude larger; see [LWS ∗
15, Fig.16].
Techniques that do not attempt to produce an approximation sufferfrom instability when computing the medial axis. Filtering tech-niques attempt to remove portions of the medial axis that do notcontribute significantly to the geometry reconstructed as the union ofmedial spheres. As shown in [MGP10, Fig.2], this works fairly wellfor inputs with little or no noise. Conversely, with large noise (and/oroutliers), these methods become mostly inappropriate; see Figure 3.
Angle filtering – θ -medial axis. One way to identify the signifi-cance of a point is the largest angle formed by the center of thecorresponding maximal sphere and two of its tangent points on theshape boundary. The θ -medial axis of [FLM03], filters out balls as-sociated with low aperture angle. While this filtering can disconnectportions of the medial axis, see Figure 3(c), the issue can be avoidedby homotopy-preserving pruning [AM97]. Distance filtering – λ -medial axis. Another metric for filteringdiscards a sphere whenever its tangent points lie on the surfacebelow a certain distance [CL05, ACK01]. As Figure 3(d) illustrates,this results in a loss of features even before all noise has beenremoved. This shortcoming makes these solutions inappropriatewhenever the input shape contains structures at different scales.
Scale filtering – σ -medial axis. The scale axis transform intro-duced by Giesen et al. [GMPW09] and extended to 3D by Mik-los et al. [MGP10] are methods built over the maximality propertyof the medial axis. First medial spheres are scaled by a given σ > / σ .This method is global as it processes the whole point set at once,and computes multiple Voronoi diagrams in the process, resulting inreduced computational efficiency – e.g. ≈ min . for a mesh with ≈ k vertices, see [MGP10, Tab.1]. The definitive version is available at wileyonlinelibrary.com ebain et al. / LSMAT: Least Squares Medial Axis Transform
Figure 4: (left) The sphere-shrinking of Ma et al. [MBC12], and(right) its resulting maximal spheres on a noisy point set.
3. Technical details
In Section 3.1, we start by reviewing the highly relevant method ofMa et al. [MBC12], which computes the medial axis by searching formaximal spheres via local analysis. We then formulate our numericaloptimization to compute medial spheres in Section 3.2.
Notation.
We are given in input an oriented point cloud { ( p n , n n ) } n = ... N drawn from the solid object O with watertightsurface/boundary ∂ O . The cloud is affected by noise with standarddeviation σ p . We give all numerical values of σ p and other parame-ters which represent distance values as percentages relative to thediagonal of the bounding box of the input shape. With M we refer tothe internal medial axis of O ; see [TDS ∗
16, Fig.2]. With s = ( c , r ) we indicate a sphere centered at c of radius r , and with t the indexof the solver iteration. The algorithm proposed by Ma et al. [MBC12] is an iterative methodthat shrinks an initially large sphere until it satisfies the medial axisproperties. Assuming p is a point on the boundary, this processleverages two properties: (cid:172) that the sphere passing through p shouldbe empty; and (cid:173) that p − c for a smooth input curve is parallel to thecontact point normal n . The sphere-shrinking algorithm is initializedwith a large sphere, passing through p , containing the entire pointset, and whose center c lies along the ray ( p , − n ) . Whenever thedistance to the closest point f ∈ ∂ O from c is less than r (i.e. thesphere is not empty), the center c is updated by computing the spheretangent to ( p , n ) and passing through f ; see Figure 4. Finally, theloop terminates when the sphere radius change across iterationsis below numerical precision. This algorithm is highly efficient aseach sample can be processed completely independently from theothers, but as it treats all points as hard constraints in a combinatorialmanner it does not cope well with noisy inputs. Our formulationborrows from this one, but generalizes it by (cid:172) reformulating it intoa continuous optimization problem, and (cid:173) making it more robust tonoise and outliers. We define a least-squares optimization capable of generating maxi-mally inscribed spheres given an oriented point cloud P . Our opti-mization energy for a sphere s consists of the combination of two terms: E medial = ω E maximal + ω E inscribed (1)We will now proceed to describe these terms in detail, and thenprovide comparison against alternative formulations in Section 4.2. Maximality.
The maximality energy creates a constant positivepressure term that tends to increase the sphere size at each iteration.By design, we define this energy to have a constant magnitude ,that is, a contribution to the optimization energy that is constantregardless of the radius of the given medial sphere. This is achievedby considering the radius at the previous optimization iteration witha constant offset ε : E maximal = || r − ( r t − + ε ) || (2) Inscription.
Consider a signed distance function (SDF) Φ ( x ) , where Φ ( x ) < x inside the shape. If the expression of this func-tion were available to us, an inscription constraint could be easilyexpressed by the inequality Φ ( c ) < − r , which we can convert inour least squares formalism as: E inscribed = R ( r + Φ ( c )) (3)where the ramp function R ( x ) = max ( x , ) only penalizes a spherewhen it violates an inscription constraint. Unfortunately, an SDFfunction for our oriented point set P is not readily available withoutcomputing a full surface reconstruction of the point cloud [BTS ∗ ∂ O – that is, as Φ →
0. For example, the MLS formulationby Kolluri [Kol08] approximates Φ as a weighted sum of locallysupported point-to-plane functions: Φ ( x ) = ∑ n ϕ n ( x , h ) n n · ( x − p n ) ∑ n ϕ n ( x , h ) (4)where in our case ϕ n ( x , h ) = ϕ ( (cid:107) x − p n (cid:107) , h ) and ϕ ( x , h ) is asmoothly decaying radial basis function [GG07] with compact sup-port [ , h ] : ϕ ( x , h ) = b (cid:16) xh (cid:17) ; b ( x ) = (cid:40) ( − x ) x < x ≥ Φ at the medial center,which could be potentially far from ∂ O , this representation is notimmediately appropriate. However, in the context of registration,several works have shown how a quadratic approximation of Φ can be built by appropriately blending point-to-point with point-to-plane distance functions [PH03, MGPG04]. In more detail, point-to-point distances are a good approximation of Φ in the far-field (i.e. for Φ (cid:29) near-field (i.e. for Φ → c n = c t − − n n ( c t − − p n ) · n n , the projection of the center on the hyperplane ofpoint p n . We can use this to interpolate between the two metrics: Φ blend ( s ) = mix ( Φ plane ( s ) , Φ point ( s ) , ϕ n ( ¯ c n , h blend )) (6)where Φ plane ( s ) = R ( r − ( p n − c ) · n n ) represents the distance fromthe sphere to a plane, Φ point ( s ) = R ( r − (cid:107) p n − c (cid:107) ) is the euclideandistance from the sphere to a point, and mix ( a , b , x ) = xa + ( − x ) b is the linear interpolation operator; see Figure 5. Based on the The definitive version is available at wileyonlinelibrary.com ebain et al. / LSMAT: Least Squares Medial Axis Transform cp p cp p c c (a) (b) (c) (d) Figure 5: (a) The function Φ halfplane uses point normals to give a signed distance value that penalizes spheres which move outside the shape,but may give nonsensical values for points belonging to parts of the surface that are not well approximated by the sphere. (b) The function Φ point does not suffer from this issue but does not enforce that the sphere should respect the orientation of the surface. (c) We mix between thetwo distances independently for each point using the distance between the point and the projection of the sphere center onto the plane definedby the point. (d) This allows us to use the sided distance only for points where it makes sense, and use the unsigned point-to-plane functionotherwise. iteration 0 iteration 4 iteration 8 iteration 16 converged Figure 6: Iterations of LSMAT optimization on an input noisy point set. We randomly initialize sphere centers and radii, and demonstrate theexcellent convergence properties of our optimization. Here the center-pin correspondence is marked by the dotted line. Notice how althoughsome spheres are initially outside, the optimization pushes them into the interior of the point set.geometry of Φ , we can then re-formulate our inscription withrespect to each oriented point p n . Analogously to Equation 4, thisenergy is accumulated over all points in P : E inscribed ≈ ∑ n ϕ ( R ( || c t − − p n || − r t − ) , h support ) (cid:124) (cid:123)(cid:122) (cid:125) (cid:54) = Φ blend ( s ) (7)The parameter h blend defines the scale used for blending betweendistance types, and h support limits how far outside of a sphere apoint may be before its contribution falls to zero. As is typicalin robust optimization (e.g. IRLS) the weights are computed withrespect to the parameters s t − at the previous iteration – inscriptionis evaluated only for points within, or in proximity of the sphere inits previous geometric configuration s t − .Figure 7: While the algorithm is initialized in a neighborhood, op-timizing for larger spheres will cause the solver to have the sphereconverge to areas of locally maximal radius. Pinned spheres.
The medial axis provides an estimate of the localthickness of the shape through its sphere local radius. However,as our variational formulation attempts to create larger spheres,nothing prevents a sphere from traveling along medial brancheswherever we have a non-vanishing medial radius gradient on M ;e.g. a sphere would slide from the tip of a cone to its base; seeFigure 7. We can avoid this issue by “pinning” medial spheres.Generalizing the exact constraints of the sphere-shrinking algorithmby [MBC12], we subject the sphere associated with a given point p to the hard constraint (cid:107) c − p (cid:107) − r ≤ d pin , which keeps the spherein contact with a sphere of radius d pin centered at p . We includethis constraint in our optimization via the penalty method [NW06],yielding a quadratic barrier energy: E pinning = R ( (cid:107) c − p (cid:107) − ( r + d pin )) (8) Optimization.
Similarly to [LCOLTE07], our optimization inducesthe definition of medial sphere as the fixed point solution of anupdate equation: s = F ( s ) = arg min s E medial + E pinning (9)Our optimization problem is quadratic, but while its energy termsare differentiable, they are not linear. Hence, we iteratively computea solution via Gauss-Newton. This requires the linearization of thearguments of the quadratic functions with respect to c and r , whichis straightforward in our setting. The definitive version is available at wileyonlinelibrary.com ebain et al. / LSMAT: Least Squares Medial Axis Transform
Implementation.
Due to the independent nature of the per-sphereoptimization, our implementation is straightforward. At each stepwe compute Jacobian matrices and residual vectors separately foreach sphere and thus are only required to solve a
DxD linear systemfor each, where D is the number of degrees of freedom of a singlesphere (3 or 4 in our experiments). This property eliminates the needfor any advanced solvers or factorizations and enables the imple-mentation to be completely parallel over the spheres, as required foreffective GPU acceleration. Additional performance optimizationcould be achieved by collecting the points for each sphere whosecontribution is non-zero using accelerating data structures such askd-trees [FBF77].
4. Results and evaluation
We show medial representations generated by our method, in both2D and 3D, throughout the paper and in Figure 16. As shown, ourmethod produces valid medial representations for a wide varietyof shapes and models, and in the presence of noise and outliers.Throughout the paper, we process all input with the same parame-ters whose values were set according to the analysis in Section 4.6.We evaluate our method in several ways. First, we consider vari-ations of the inscription and maximality energy, and show howalternative formulations fail. Second, we evaluate the performanceof our algorithm against ground truth on synthetic benchmarks con-taminated by noise and outliers. Third, we compare our generatedmedial representations against the state-of-the-art in both 2D and3D, again in the presence of noise. Finally, we provide an analysisof our algorithm parameters.
We consider two modified formulations of Equation 7 in which weeither penalize the squared distance to points, or the squared distanceto half-planes as opposed to our blended formulation. We investigatethis behavior by randomly initializing the algorithm, and snapping ϕ n ( ¯ c n ) to either zero or one for all points. When ϕ n ( ¯ c n ) =
1, point-to-point energy is used and the algorithm attempts to make spheresempty in a least square sense. However, nothing prevents the opti-mization from generating maximal spheres outside the shape in theambient space. When ϕ n ( ¯ c n ) =
1, point-to-plane energy is used, and Φ halfplane Φ point LSMAT Figure 8: Qualitative evaluation of inscription energy variants. Figure 9: Qualitative evaluation of maximality energy variants.the algorithm attempts to make half-spaces empty . As illustratedin Figure 8 and the supplemental video , this results in difficultiesfor the algorithm in dealing with sharp concavities. In the exampleabove, we expect all centers to cluster to one of the two centers,but spheres in the neighborhood of the sharp concavity might readthe normal of a point on the opposite side, with a halfplane request-ing the radius of a sphere intersecting with it to be significantlysmaller. This problem is caused by the fact that point-to-point andpoint-to-plane energies approximate the squared SDF of a pointset respectively in the far and near field. Our LSMAT formulationrespects this geometric property, and deals with both issues at once.
We consider two alternative formulations of the maximality con-dition in our optimization: penalizing the squared inverse of thesphere radius, expressed as (cid:107) / r (cid:107) , or directly specifying R max ,the maximum size of a medial sphere, and optimizing (cid:107) r − R max (cid:107) .While intuitively these are potentially feasible solutions, they incura significant limitation: the amount of “pressure” a medial spherewill apply will be dependent on its size. That is, small spheres willbe associated with a higher energy level. For example, using the firstformulation, as r → −∞ ; as illustratedin Figure 9, this can cause spheres in the neighborhood of smallfeatures to bulge out. Our solution, detailed in Equation 2, does notencounter this problem, as our energy is constant regardless of thesphere size.Figure 10: As the optimization of Equation 9 is executed, the av-erage/max ground truth errors converge to the maximum precision.This plot illustrates the error for the iterations visualized in Figure 6. The definitive version is available at wileyonlinelibrary.com ebain et al. / LSMAT: Least Squares Medial Axis Transform
LSMAT %% Figure 11: (top) Quantitative evaluation on state-of-the-art methodsshowing increasing error as we increase the noise level. (bottom)Qualitative results corresponding to the two dashed lines in the plot.
We consider ground truth geometry polluted by noise, and evaluatethe quality of an extracted axis by computing the distance of eachmedial center to the closest point on the ground truth axis ˜ M : E avg = N ∑ n arg min ˜ c n ∈ ˜ M (cid:107) c n − ˜ c n (cid:107) (10) E max = max n arg min ˜ c n ∈ ˜ M (cid:107) c n − ˜ c n (cid:107) (11)We compute ˜ M via the medial axis module of the python skimage package, derived from the ridges of the distance transform from theground truth boundary ∂ O . As we discretize images with a boundingbox at a 1024 resolution, the “numerical precision” of the metricsabove is of one pixel. To obtain scale invariance, we report theseerrors in relation to the diagonal of the bounding box. In Figure 10,we visualize the convergence of our iterative optimization scheme. Through the metric from Section 4.3, we quantitatively evaluate theperformance of LSMAT against local filtering methods, as well asthe global scale axis transform. Figure 3 shows how neither distancenor angle filtering is very effective; hence we employ a compoundvariant where we first filter by distance with λ > σ p , and then byangle θ > ◦ . For the scale axis, we set σ = .
3, and for ourmethod, we use our default parameters. We gradually increment thelevel of noise, and plot the corresponding error metric in Figure 11,as well as a few example frames from the plot. Note that LSMATcorrectly captures the shape of the ground truth boundary, whereasall other local methods fail.
We qualitatively evaluate the performance of LSMAT in threedimensions by comparing our results to those generated by the SAT [MGP10], QMAT [LWS ∗ σ p ∈ [ , ] relative to the diagonal boundingbox. On the fertility model, LSMAT produces a convincing medialaxis representation even when the input oriented point cloud is af-fected by extreme noise – i.e. with a magnitude close to the one ofthe local feature size . For the vase model, LSMAT still producesa smoother and more faithful medial axis representation than theexisting state-of-the-art. Even in the presence of minor amounts ofnoise (second row), LSMAT faithfully produces a medial axis repre-sentation with a smooth surface, and whose skeleton resembles thatof the uncorrupted model. The results are particularly encouragingon the horse dataset, where we attempted to use a very small numberof target primitives for QMAT ( ≈ ≈ . Our algorithm depends on five parameters: the kernel sizes h blend and h support , the relative weight ω / ω , pinning distance d pin , andradius expansion constant ε . The effect of varying the first fourgiven different levels of noise is shown in Figure 13. Note we onlyconsider the ratio between the ω ∗ , as the two energies balance eachother. As h blend increases, spheres near sharp concave corners beginto shrink as they eventually use the point-to-plane distance for allpoints in their neighborhood. When h support grows significantly be-yond the local feature size, the MLS formulation is no longer able torecover a meaningful surface. As ω / ω increases and the relativeimportance of maximality versus inscription increases, spheres ex-pand until they no longer form a faithful medial representation andultimately escape the point cloud. Finally, as d pin grows, the spheresare allowed to slide further from their starting positions towards ar-eas of locally maximal radius. We assume that estimates of the inputnoise characteristics are available, and choose our default value forModel ( σ p ) LSMAT Sphere Shrinking
10k Spheres Time Iterations Time IterationsFertility (0%) 11.8s 40 0.30s 9Fertility (1%) 11.6s 40 0.35s 11Fertility (2%) 12.0s 40 0.33s 10Fertility (5%) 16.2s 40 0.34s 10Vase (0%) 11.2s 70 0.16s 9Vase (1%) 12.1s 70 0.17s 10Vase (2%) 11.8s 70 0.16s 9Vase (5%) 6.1s 70 0.16s 9Horse (2%) 6.8s 60 0.16s 9Table 1: Run-times for LSMAT and Sphere Shrinking results pre-sented in Figure 12. Both algorithms are GPU accelerated and runon the same hardware. Note that our efficiency claims are madewith respect to QMAT and the scale axis transform, for which GPUacceleration is not availible.
The definitive version is available at wileyonlinelibrary.com ebain et al. / LSMAT: Least Squares Medial Axis Transform
Input QMAT Scale Axis Sphere Shrinking LSMAT F e r tilit y V a s e H o r s e Figure 12: Qualitative evaluation on state-of-the-art methods as we increase the noise level. The QMAT [LWS ∗
15] and SAT [MGP10] in thefirst two columns are global methods that assume a watertight surface, while the sphere-shrinking [MBC12] and the LSMAT proposed hereare local methods.each of these four parameters based on empirically derived linearfunctions of σ p ; see Appendix B for details. Through experimentswe observed that it is sufficient to choose a constant value for thefifth parameter ε , as different values of it merely change the optimalrelation between ω / ω and σ p . For all our experiments we use ε = The definitive version is available at wileyonlinelibrary.com ebain et al. / LSMAT: Least Squares Medial Axis Transform h blend h support ω / ω d pin σ p = % σ p = % σ p = % σ p = % σ p = % σ p = % σ p = % σ p = % σ p = % σ p = % σ p = % σ p = % Figure 13: A qualitative analysis of parameters in our algorithm. Each quadrant shows the result of sweeping a parameter (horizontal) fordifferent noise values (vertical). The highlighted images represent the "default" parameter choice as defined in Section 4.6. The right columnin each quadrant shows an unreasonably high choice to demonstrate that there is an upper bound for each parameter.
5. Future works
Expressing the Medial Axis Transform as a non-linear least squaresproblem opens up several interesting avenues for future research.
Robustness to outliers – Figure 14.
Least squares problems canbe interpreted as a maximum-likelihood (ML) estimation given a
Gaussian probability distribution of the noise variables. If the inputis corrupted by other forms of noise, one could replace
Gaussian with other error distributions, and derive the corresponding MLoptimization scheme. For example, if we assume the probabilitydistribution of the noise to be Laplacian, the least-squares problemswould simply be transformed into an (cid:96) (i.e. least norm) optimiza-tion. However, these type of problems can still be computed withGauss-Newton type methods by using iteratively re-weighted leastsquares (IRLS) techniques [DDFG10]. As illustrated in Figure 14this results in an IR-LSMAT algorithm that can cope with significant amounts of outliers. While these results are promising, the conver-gence speed of the optimization is severely reduced, as a muchsmaller value for ε was needed to produce these results. The gen-eralizability of LSMAT to (cid:96) p robust norms [BTP13] is particularlyinteresting. More specifically, consider the optimization problemconsisting only of the following energy:arg min c , r ∑ n |(cid:107) c − p n (cid:107) − r | p . (12)For p = p → symmetry set , a superset of the medial axis;see Figure 2(d) and [TDS ∗
16, Sec. 2.1.4]. How to exploit (12) to
The definitive version is available at wileyonlinelibrary.com ebain et al. / LSMAT: Least Squares Medial Axis Transform
Figure 14: Iteratively Reweighted LSMAT and its ability to cope with an increasing number of outliers (as % of input point set).efficiently compute the MAT of a point set is an interesting venuefor future works.
Smoothness priors – Figure 15.
Our pinning formulation is highlyefficient, but its local nature can be also considered a intrinsic limi-tation. For example, in the noisy maple leaf example in Figure 16,at times the estimated medial spheres could be over and/or under-estimated in size, resulting in medial centers that do not necessarilysample the underlying piecewise-smooth manifold of the MAT. How-ever, if our input is a sampling of a smooth surface ∂ O , then weknow that [SP08]: (cid:172) the medial centers c ∗ should be lying on apiecewise smooth manifold, and (cid:173) the sphere radius function onthis manifold should vary smoothly. Another desirable character-istic might be to have a uniform sampling of this manifold. Wecan convert these priors into least-squares energies, resulting in a initialization iteration Figure 15: Optimizing LSMAT centers placement. To highlight thesmoothness of the resulting shape, we only display the input pointset overlaid to the initialization. maximum a-posteriori optimization. In Figure 15 we illustrate afew iterations of this optimization on the moon shape, where the“pinning” constraints from Section 3.2 have been disabled. While theoptimization behaves as expected, this suffers similar shortcomingsto those illustrated in Figure 7, and would result in a single sphere ifexecuted for t → ∞ . Modifying the variational LSMAT formulationto obtain a regularly sampled distribution of medial centers is aninteresting venue for future works. Optimization acceleration.
In our experiments, we initialize theoptimization with random sphere positions and radii. Nonetheless,given how a smooth object is composed of smooth piecewise man-ifolds M , and a smooth radius function R defined thereon, onecould easily envision a locally bootstrapped version of the algo-rithm, where unsolved medial balls are initialized with the c ∗ , r ∗ of their neighbors – which could also be re-interpreted as a multi-scale solver. Notice that techniques such as the sphere-shrinkingalgorithm from [MBC12] do not permit this type of acceleration. Shape approximation vs. reconstruction.
A number of methodsleverage the medial axis for interpolatory reconstruction [Dey06],and our work is a stepping stone towards the creation of approximat-ing reconstruction [BTS ∗
16] algorithms based on the medial axis.Methods such as Medial Meshes [SCYW16] and QMAT [LWS ∗ already avail-able, and attempt to extract its approximation via swept-spheres.Conversely, our work could be extended to directly compute a reconstruction of the input point cloud by minimizing data fit-ting metrics based on Hausdorff distances [SCYW16] or sphericalquadrics [LWS ∗ Orienting a point set.
Finally, while our formulation is based onan oriented point set and an approximation of its SDF in the near/farfield, an interesting variant of our algorithm could consider an un-oriented point set, where the quantities to be optimized for wouldbe the radii [ r , r ] , and contact plane [ k , n ] of a pair of twin spheres .The contact point t should then be then be optimized on the manifold ∂ O , while the complementary outside/inside label of each spherebe optimized to result in a smooth signing of the environment space.This approach would then provide a “medial axis” analogous torecent efforts in variational reconstruction of non-oriented pointclouds [MDGD ∗ The definitive version is available at wileyonlinelibrary.com ebain et al. / LSMAT: Least Squares Medial Axis Transform
6. Conclusions
We have introduced the Least Squares Medial Axis Transform, orLSMAT, a continuous relaxation of the medial axis transform that isnot only stable, but also robust to high levels of noise even thoughit is based solely on local optimization. While in most of the paperwe visualized the generated maximal spheres covering more orless the entire shape, we would like to remind the reader that thealgorithm operates on each sphere independently ; the algorithm istherefore trivially parallelizable and particularly suitable for GPUimplementations. Our method produces results on noisy inputs thatstate-of-the-art methods fail to handle, without a reliance on ad-hoc postprocessing. Our approach is efficient, parallelizable, andtherefore suitable for real time applications where reliable medialrepresentations are required, and where captured inputs are likely tobe noisy.
Acknowledgements
We would like to express our gratitude to Mathieu Desbrun,Leonidas Guibas, Justin Solomon, Keenan Crane, and especiallySofien Bouaziz for their ideas and suggestions. In our experimentswe use models provided courtesy of [CGF09] and the AIM@SHAPEShape Repository.
References [AB99] A
MENTA
N., B
ERN
M.: Surface reconstruction by voronoi filter-ing.
Discrete & Computational Geometry (1999). 3[ACK01] A
MENTA
N., C
HOI
S., K
OLLURI
R. K.: The power crust. In
Proceedings of the ACM Symposium on Solid Modeling and Applications (Ann Arbor, Michigan, USA, 2001). 3[AM97] A
TTALI
D., M
ONTANVERT
A.: Computing and simplifying 2dand 3d continuous skeletons.
Computer vision and image understanding (1997). 3[BA92] B
RANDT
J. W., A
LGAZI
V.: Continuous skeleton computationby voronoi diagram.
CVGIP: Image Understanding (1992). 2[Blu67] B
LUM
H.:
A transformation for extracting new descriptors ofshape . M.I.T. Press, 1967. 1[BP07] B
ARAN
I., P
OPOVI ´C
J.: Automatic rigging and animation of 3dcharacters.
ACM Trans. on Graphics (Proc. of SIGGRAPH) (2007). 1[BTP13] B
OUAZIZ
S., T
AGLIASACCHI
A., P
AULY
M.: Sparse iterativeclosest point.
Computer Graphics Forum (Proc. of SGP) (2013). 9[BTS ∗
16] B
ERGER
M., T
AGLIASACCHI
A., S
EVERSKY
L., A
LLIEZ
P.,G
UENNEBAUD
G., L
EVINE
J., S
HARF
A., S
ILVA
C.: A survey of surfacereconstruction from point clouds.
Proc. Eurographics (State of the ArtReports) (2016). 4, 10[CGF09] C
HEN
X., G
OLOVINSKIY
A., F
UNKHOUSER
T.: A benchmarkfor 3d mesh segmentation. In
Acm transactions on graphics (tog) (NewOrleans, Louisiana, USA, 2009). 11[CL05] C
HAZAL
F., L
IEUTIER
A.: The λ -medial axis. Graphical Models (2005). 3[DDFG10] D
AUBECHIES
I., D E V ORE
R., F
ORNASIER
M., G
ÜNTÜRK
C. S.: Iteratively reweighted least squares minimization for sparse recov-ery.
Communications on pure and applied mathematics (2010). 9[Dey06] D EY T. K.:
Curve and surface reconstruction: algorithms withmathematical analysis . Cambridge University Press, 2006. 10[FBF77] F
RIEDMAN
J. H., B
ENTLEY
J. L., F
INKEL
R. A.: An algorithmfor finding best matches in logarithmic expected time.
ACM Transactionson Mathematical Software (TOMS) (1977). 6 [FLM03] F
OSKEY
M., L IN M. C., M
ANOCHA
D.: Efficient computationof a simplified medial axis. In
Proceedings of the ACM Symposium onSolid Modeling and Applications (Seattle, Washington, USA, 2003). 3[GG07] G
UENNEBAUD
G., G
ROSS
M.: Algebraic point set surfaces.
ACM Transactions on Graphics (TOG) (2007). 4[GH97] G
ARLAND
M., H
ECKBERT
P. S.: Surface simplification usingquadric error metrics. In
Proc. of ACM SIGGRAPH (Los Angeles, Cali-fornia, USA, 1997). 3[GMPW09] G
IESEN
J., M
IKLOS
B., P
AULY
M., W
ORMSER
C.: Thescale axis transform. In
Proceedings of the Twenty-fifth Annual Sympo-sium on Computational Geometry (Aarhus, Denmark, 2009). 3[JKT13] J
ALBA
A. C., K
USTRA
J., T
ELEA
A. C.: Surface and curveskeletonization of large 3d models on the gpu.
IEEE Transactions onPattern Analysis and Machine Intelligence (2013). 3[Kol08] K
OLLURI
R.: Provably good moving least squares.
ACM Trans-actions on Algorithms (TALG) (2008). 4[LCOLTE07] L
IPMAN
Y., C
OHEN -O R D., L
EVIN
D., T AL -E ZER
H.:Parameterization-free projection for geometry reconstruction.
ACM Trans.on Graphics (Proc. of SIGGRAPH) (2007). 5[LWS ∗
15] L I P., W
ANG
B., S UN F., G UO X., Z
HANG
C., W
ANG
W.: Q-mat: computing medial axis transform by quadratic error minimization.
ACM Trans. Graph. (2015). 3, 7, 8, 10[MBC12] M A J., B AE S. W., C
HOI
S.: 3d medial axis point approxima-tion using nearest neighbors and the normal field.
The Visual Computer (2012). 2, 3, 4, 5, 7, 8, 10[MDGD ∗
10] M
ULLEN
P., D E G OES
F., D
ESBRUN
M., C
OHEN -S TEINER
D., A
LLIEZ
P.: Signing the unsigned: Robust surface reconstruction fromraw pointsets.
Computer Graphics Forum (Proc. of SGP) (2010). 10[MGP10] M
IKLOS
B., G
IESEN
J., P
AULY
M.: Discrete scale axis repre-sentations for 3d geometry.
ACM Trans. Graph. (2010). 3, 7, 8[MGPG04] M
ITRA
N. J., G
ELFAND
N., P
OTTMANN
H., G
UIBAS
L.:Registration of point cloud data from a geometric optimization perspective.In
Proc. of SGP (Nice, France, 2004). 4[MHR ∗
16] M
USIALSKI
P., H
AFNER
C., R
IST
F., B
IRSAK
M., W
IMMER
M., K
OBBELT
L.: Non-linear shape optimization using local subspaceprojections.
ACM Trans. on Graphics (Proc. of SIGGRAPH) (2016). 1[NW06] N
OCEDAL
J., W
RIGHT
S. J.:
Sequential quadratic programming .Springer, 2006. 5[PH03] P
OTTMANN
H., H
OFER
M.: Geometry of the squared distancefunction to curves and surfaces. In
Visualization and mathematics III .Springer, 2003, pp. 221–242. 4[SBdB16] S
AHA
P. K., B
ORGEFORS
G., DI B AJA
G. S.: A survey onskeletonization algorithms and their applications.
Pattern RecognitionLetters (2016). 2[SCYW16] S UN F., C
HOI
Y. K., Y U Y., W
ANG
W.: Medial meshes:a compact and accurate representation of medial axis transform.
IEEETransactions on Visualization and Computer Graphics (2016). 3, 10[SP08] S
IDDIQI
K., P
IZER
S.:
Medial representations: mathematics,algorithms and applications . Springer, 2008. 10[SSCO08] S
HAPIRA
L., S
HAMIR
A., C
OHEN -O R D.: Consistent meshpartitioning and skeletonisation using the shape diameter function.
TheVisual Computer (2008). 1, 3[TD17] T
SOGKAS
S., D
ICKINSON
S.: Amat: Medial axis transform fornatural images. In (Venice, Italy, 2017). 1[TDS ∗
16] T
AGLIASACCHI
A., D
ELAME
T., S
PAGNUOLO
M., A
MENTA
N., T
ELEA
A.: 3d skeletons: a state-of-the-art report.
Proc. Eurograph-ics (State of the Art Reports) (2016). 2, 4, 9[TGB13] T
HIERY
J.-M., G UY E., B
OUBEKEUR
T.: Sphere-meshes:Shape approximation using spherical quadric error metrics.
ACM Trans.Graph. (2013). 3
The definitive version is available at wileyonlinelibrary.com ebain et al. / LSMAT: Least Squares Medial Axis Transform [TGBE16] T
HIERY
J.-M., G UY E., B
OUBEKEUR
T., E
ISEMANN
E.:Animated mesh approximation with sphere-meshes.
ACM Trans. Graph. (2016). 3[TPT16] T
KACH
A., P
AULY
M., T
AGLIASACCHI
A.: Sphere-meshes forreal-time hand modeling and tracking.
ACM Transactions on Graphics(Proceedings of SIGGRAPH Asia) (2016). 1
Appendix A:
Gradients of Energy ComponentsGradients are written in the form ∇ s f ( s ) = (cid:104) ∇ c f ( s ) , ∂ f ( s ) ∂ r (cid:105) . Inscription Term.
Given oriented surface point ( p n , n n ) : ∇ s Φ plane ( s ) = H ( Φ plane ( s ))[ − n n , ] (13) ∇ s Φ point ( s ) = H ( Φ point ( s )) (cid:20) p n − c (cid:107) p n − c (cid:107) , (cid:21) (14) Maximality Term. R maximal = r − ( r t − + ε ) (15) ∇ s R maximal = [ , ] (16) Pinning Term.
Given pin point p : R pinning = R ( (cid:107) c − p (cid:107) − ( r + d pin )) (17) ∇ s R pinning = H ( R pinning ) (cid:20) c − p (cid:107) c − p (cid:107) , − (cid:21) (18) Appendix B:
Empirical Parameter DefaultsDefault value of ω / ω : ω / ω = . σ p + .
02 (19)Default value of h blend and h support (both the same): h blend = . σ p + .
49 (20)Default value of d pin : d pin = . σ p (21) The definitive version is available at wileyonlinelibrary.com ebain et al. / LSMAT: Least Squares Medial Axis Transform
Figure 16: A gallery of LSMAT results with varying levels of noise on shapes with complex topology and varying feature size. In the calloutfor the octopus, notice how the MLS kernel overlaps nearby surfaces, yet the algorithm can cope by producing erroneously located medialspheres with zero radius.
The definitive version is available at wileyonlinelibrary.com ebain et al. / LSMAT: Least Squares Medial Axis Transform σ p = 0% σ p = 0.1% σ p = 1% σ p = 2% k = 10 -6 k = 10 -5 k = 10 -4 Figure 17: Parameter sweep for QMAT. The X axis represents noise, while the Y axis shows results for different values of the parameter k . Wefound that k had little effect on the noise tolerance of the algorithm. For the left two columns, the noise values are the minimum and maximumnoise levels shown in the original publication, which yield good results. However, for the higher noise values that we test against QMAT failsto produce a useful medial representation. As noted in Section 2.1, this is a known and acknowledged limitation for this family of methods. σ p = 0% σ p = 0.1% σ p = 1% σ p = 2% σ = 1.1 σ = 1.3 σ = 1.5 Figure 18: Parameter sweep for Scale Axis. The X axis shows the same noise levels used in Figure 17, while the Y axis shows results fordifferent values of the scale parameter σ . As expected, this method performs well in areas where the local feature size is much larger than thenoise. We show the intermediate up-scaled spheres to demonstrate this more clearly. The Scale Axis Transform is in some cases able to recovera useful medial representation even when the unfiltered MA is highly corrupted. The missing images are due to the implementation failing tocomplete within the allowed time for that combination of parameters.. As expected, this method performs well in areas where the local feature size is much larger than thenoise. We show the intermediate up-scaled spheres to demonstrate this more clearly. The Scale Axis Transform is in some cases able to recovera useful medial representation even when the unfiltered MA is highly corrupted. The missing images are due to the implementation failing tocomplete within the allowed time for that combination of parameters.