Nonlinear Spectral Geometry Processing via the TV Transform
NNonlinear Spectral Geometry Processing via the TV Transform
MARCO FUMERO,
Sapienza University of Rome
MICHAEL MÖLLER,
Siegen University
EMANUELE RODOLÀ,
Sapienza University of Rome
We introduce a novel computational framework for digital geometry pro-cessing, based upon the derivation of a nonlinear operator associated tothe total variation functional. Such an operator admits a generalized notionof spectral decomposition, yielding a convenient multiscale representationakin to Laplacian-based methods, while at the same time avoiding unde-sirable over-smoothing effects typical of such techniques. Our approachentails accurate, detail-preserving decomposition and manipulation of 3Dshape geometry while taking an especially intuitive form: non-local seman-tic details are well separated into different bands, which can then be filteredand re-synthesized with a straightforward linear step. Our computationalframework is flexible, can be applied to a variety of signals, and is easilyadapted to different geometry representations, including triangle meshesand point clouds. We showcase our method through multiple applications ingraphics, ranging from surface and signal denoising to enhancement, detailtransfer, and cubic stylization.CCS Concepts: •
Computing methodologies → Shape analysis .Additional Key Words and Phrases: total variation, spectral geometry
ACM Reference Format:
Marco Fumero, Michael Möller, and Emanuele Rodolà. 2020. Nonlinear Spec-tral Geometry Processing via the TV Transform.
ACM Trans. Graph.
39, 6, Ar-ticle 199 (December 2020), 16 pages. https://doi.org/10.1145/3414685.3417849
In the last decades, the computer graphics community has witnesseda thriving trend of successful research in computational spectralgeometry. Spectral approaches have met with considerable tractiondue to their generality, compact representation, invariance to datatransformations, and their natural interpretation relating to classicalFourier analysis. Despite these benefits, the vast majority of currentspectral methods rely upon the construction of smooth basis func-tions (obtained, for instance, as the eigenfunctions of the Laplaceoperator), leading to significant detail loss whenever the signal to berepresented is not smooth. When this signal represents geometricinformation (e.g., the ( x , y , z ) coordinates of mesh vertices), this in-evitably leads to a “lossy” geometric encoding, typically manifest asover-smoothing, vertex collapse, or spurious oscillations (so calledGibbs phenomenon) not appearing in the original geometry. Authors’ addresses: Marco Fumero, [email protected], Sapienza University ofRome, Michael Möller, [email protected], Siegen University, EmanueleRodolà, [email protected], Sapienza University of Rome,Permission to make digital or hard copies of all or part of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage and that copies bear this notice and the full citationon the first page. Copyrights for components of this work owned by others than ACMmust be honored. Abstracting with credit is permitted. To copy otherwise, or republish,to post on servers or to redistribute to lists, requires prior specific permission and/or afee. Request permissions from [email protected].© 2020 Association for Computing Machinery.0730-0301/2020/12-ART199 $15.00https://doi.org/10.1145/3414685.3417849
Fig. 1. Our method allows us to process the geometric details of a givenshape into well separated spectral bands having different semantics. Eachband corresponds to a different geometric feature, possibly not localized inspace. In this example, the star decorations are removed from the knot via asimple low-pass filter in the TV spectrum. The method is rotation-invariant,as it does not depend on the orientation of the knot in 3D space.
Such artifacts derive from the motivations underlying the seminalwork of Taubin [1995] and follow-up [Desbrun et al. 1999], whichwere posed as an extension of the diffusion techniques from imageprocessing to mesh smoothing applications. Today, the key driver be-hind their continued adoption lies in the convenient representation,which enabled significant gaps in performance in many relevantproblems (e.g., shape correspondence [Melzi et al. 2019; Ovsjanikovet al. 2012], vector field processing [Brandt et al. 2017] among manyothers), while still suffering from often undesirable effects.The present work finds its motivation in the fact that detail preser-vation is a strict necessity in a wide class of applications in graphics.We start from the observation that, similar to most natural images,shape geometry often has a sparse gradient. For a scalar function f ,this quantity can be measured via the total variation functional: TV ( f ) = ∫ ||∇ f ( x )|| dx , (1)which, for coordinate functions, captures the amount of sharp geo-metric variations in a given object. Leveraging recent progress inimage processing [Burger et al. 2016; Gilboa 2013], we use the abovefunctional to derive a nonlinear operator, the subdifferential ∂ TV ,which has numerous piecewise-constant functions as eigenfunctions .Using the TV functional, we provide a spectral decomposition( analysis ) algorithm for functions f on a manifold which, similarlyto Fourier decompositions, carries a canonical multiscale ordering,providing a sharp separation of the source signal into different bands;as a direct consequence to the sparsity of the gradient, frequencybands in the spectral domain assume a semantic connotation, andcan be filtered and processed following classical signal processingparadigms (see Figure 1 for an example). Importantly, since TV ( f ) ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020. a r X i v : . [ c s . G R ] S e p remains finite for a wide variety of discontinuous f (with an ap-propriate adaptation of the definition (1) to be discussed in Eq. (4)),discontinuities are preserved by the transform, and is therefore well-suited for representing jumps and sharp high-frequency variations.The inverse transform ( synthesis ) is realized by a single linear step.Our approach shares some common traits with Laplacian meshprocessing, for which it provides a nonlinear alternative; in particu-lar, we show that both settings can be derived as specific instancesof a more general model. We demonstrate the versatility of ourframework by applying it with success to several applications ingeometry processing and graphics, relating its multiple incarna-tions to existing ad-hoc approaches, and illustrating its performanceacross different settings and shape representations. With this paper we introduce a new framework for geometry pro-cessing that weaves discontinuity-preserving regularization intothe fabric of a spectral approach. Our key contributions are: • We introduce a general spectral framework to analyze andprocess non-smooth signals defined on surfaces, as well as thesurfaces themselves, which fully preserves geometric detail; • We adopt a fast, unconditionally stable algorithm to solve theresulting nonlinear problem on non-Euclidean domains, withlimited dependence on a few parameters; • We analyze the evolution of surfaces along the flow inducedby the operator, and relate it to existing geometric flows; • We explore a wide range of possible applications in shapeanalysis and modelling, where we compare to state-of-the-artmethods in the respective domains, demonstrating productionquality results in many cases.
From a general perspective, ours is a filtering approach that incorpo-rates an edge-preserving diffusion process within a spectral decom-position framework, and draws inspiration from recent progress innonlinear image processing. In the following, we cover the worksthat more closely relate with these themes.
Our methodology falls within the class of energy minimizationmethods for geometry processing, and carries a natural dual inter-pretation as a diffusion-based approach.
Diffusion-based methods.
Following the work of Taubin [1995],classical editing tasks such as surface smoothing, denoising andenhancement have been traditionally phrased as variants of thediffusion PDE ∂ t u = ∆ M u . Here, u is a signal defined on the sur-face M , and ∆ M is its Laplace-Beltrami operator. For example,feature-aware denoising driven by local curvature has been ad-dressed by modifying the metric on M (and hence its Laplacian),leading to a rich production of anisotropic techniques [Bajaj and Xu2003; Clarenz et al. 2000; Desbrun et al. 2000; Tasdizen et al. 2002].Higher-order filtering with volume preservation to prevent surfaceshrinking was proposed in [Desbrun et al. 1999]; the same workclarified the well known duality relating the diffusion to mean cur-vature flow, and showed how the latter notion leads to a more robust formulation for denoising applications. This was later modified in[Kazhdan et al. 2012] to resolve instabilities of the flow.Our approach also involves a diffusion process, and comes withan associated geometric flow. We identify two key differences overprior work: (1) in our case, the operator involved is nonlinear, and(2) diffusion takes place within the level sets of the signal u ratherthan across the boundary of the level sets. As a consequence, discon-tinuities are not smeared throughout the diffusion process, whilecontrast is reduced. In the inset figure, the source signal (left) isprocessed via the classical “horizon-tal” isotropic diffusion (middle) andwith our “vertical” diffusion (right);in the latter case, the width remains constant while the amplitudedecreases. When the signal encodes surface geometry, this impliesthat sharp features are preserved. Energy minimization methods.
Other approaches are based onformulating editing operations as the minimization of an energy.Although mathematically any variational problem can be rephrasedas a diffusion problem, choosing the former enables the implemen-tation of a richer class of algorithms. Such methods lift ideas fromimage processing to the surface domain. For example, in [He andSchaefer 2013] mesh denoising is done by minimizing the L gradi-ent norm of the vertex positions, thus suppressing low-amplitudedetails, as done in [Xu et al. 2011] for images; extensions to meshesof bilateral [Tomasi and Manduchi 1998] and mean shift [Comani-ciu and Meer 2002] filtering were proposed in [Fleishman et al.2003; Solomon et al. 2014]; a variant of [Zhang et al. 2014] wasapplied to the normal vector field in [Wang et al. 2015]; shock fil-ters [Osher and Rudin 1990] were applied to surface geometry forfeature sharpening and enhancement in [Prada and Kazhdan 2015].Gradient-domain mesh processing [Chuang et al. 2016] is also arecent direction inspired by success in the image domain, allowingexplicit prescription of target gradient fields.The above testifies to a wealthy literature of nonlinear edge pre-serving strategies for mesh filtering. The main drawback of thesemethods lies in their lack of generality. Here we provide a generalframework to perform multiscale analysis and synthesis of arbitrarysignals defined on geometric domains (meshes and point cloudsalike); while denoising, fairing and enhancement can be easily real-ized within this framework, we are not limited to these applicationsand provide a more flexible setting for geometry processing. TV-based methods.
TV and the strictly related Mumford-Shahfunctional have been considered before for different tasks in geom-etry processing, e.g. surface reconstruction [Liu et al. 2017]. Zhanget al. [2015] consider an extension of the ROF model [Rudin et al.1992] for denoising surfaces, which they solve with an augmentedLagrangian method followed by a geometry reconstruction step[Sun et al. 2007]. Zhong et al. [2018] follow a similar approach, withthe inclusion of an additional regularizer defined as the Laplacianenergy of the normal coordinates, which helps dampen the staircaseeffect that is observed in [Zhang et al. 2015]. These methods arerelated to ours, as they involve the minimization of the TV of thesurface normal field; we also consider this kind of energy, but inaddition, we study the geometric flow associated to it. Further, ourapproach is not restricted to normal fields, but allows processing
ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020. onlinear Spectral Geometry Processing via the TV Transform • 199:3 arbitrary surface signals. Tong and Tai [2016] consider the Mumford-Shah functional as a measure of sharp geometric variation for agiven signal. Their focus is on the detection of features lines (i.e.the boundaries between neighboring segments) for scalar functionssuch as mean curvature and color texture, rather than geometryprocessing. The MS functional has also been considered in [Bonneelet al. 2018] in its Ambrosio-Tortorelli approximation. The authorsextend the range of possible applications of [Tong and Tai 2016] byoperating on the surface normal field, showing results in inpainting,denoising and segmentation tasks. In addition to lacking an analysisof the TV normal flow, none of these methods provide a spectral framework based on the TV functional for geometry processing.
As shown in Figure 1, our multiscale approach provides a sharpseparation of an input signal into “feature bands”. We achieve thisby means of a generalized spectral decomposition. Spectral geom-etry processing was introduced by Taubin [1995] and further de-veloped in [Lévy 2006; Rong et al. 2008; Sorkine 2005; Zhang et al.2010] for tasks of shape reconstruction, modelling, and deformationtransfer among others. These works generalize Fourier analysis tosurfaces by observing thatthe Laplacian eigenvectorsform an orthogonal set ofsmooth basis functions witha canonical ordering (low tohigh frequencies). This basisis optimal for representing smooth functions compactly [Aflalo et al.2014], but it gives rise to artifacts in the presence of sharp variations,thus reducing the benefits of working with a spectral representation;in the inset, we show the typical blurring effect (middle) and Gibbsoscillations (right) due to the jumps in the source signal (left), nothandled well by the smooth Laplacian eigenbasis.Variants of the Laplace operator add flexibility and feature aware-ness, including anisotropic [Andreux et al. 2015] and extrinsic [Liuet al. 2017] variants, see [Wang and Solomon 2019] for a recentsurvey. Sparsity has been analyzed in [Neumann et al. 2014], whereLaplacian eigenfunctions are modified to have limited support bypenalizing their L norm. We differ in that we consider signals withsparse gradients instead. All these operators are linear, and theirassociated energies promote smoothness. Our operator is nonlin-ear, hence more powerful, and its eigenfunctions are well suited torepresent jump discontinuities in the signal. This yields a cleanerseparation of different features into different bands, in contrast toLaplacian-based representations, which spread individual featuresacross the entire spectrum. Finally, similar to prior research, our approach is also based on ideasoriginating from image processing. In this field, nonlinear varia-tional methods have replaced linear methods due to their higherprecision in approximating natural phenomena. To date, a majorrole has been played by the total variation (TV) energy and its min-imization, originally exploited by Rudin et al. [1992] with the ROFmodel for image denoising. Other uses of the TV functional span image deblurring, inpainting, interpolation, image decomposition,super-resolution and stereovision, to name but a few; for a generalintroduction to the topic see [Chambolle et al. 2010]. Recently, aspectral theory related to the gradient flow of this functional (c.f.Bellettini et al. [2002]) on eigenfunctions, was developed by Gilboa[2013; 2014] and generalized in [Burger et al. 2016; Gilboa et al.2015] with applications to image processing, allowing for the analy-sis of images at different scales by exploiting the edge-preservingproperty of the total variation. Here we propose to lift this latterframework from the Euclidean domain to surfaces, put it in relationto existing geometry processing models, and demonstrate severalpossible applications in the area of 3D graphics.
We model shapes as compact and connected Riemannian mani-folds M embedded in R , possibly with boundary ∂ M , and withtangent bundle T M = (cid:208) p ∈M T p M . The manifold is equippedwith a metric, i.e., an inner product defined locally at each point, ⟨· , ·⟩ p : T p M × T p M → R , making it possible to compute lengthsand integrals on M . The intrinsic gradient ∇ M , defined in terms ofthe metric, generalizes the notion of gradient to manifolds; similarly,the intrinsic divergence div M can be defined as the negative adjointof the gradient, assuming Neumann boundary conditions: ∫ M div M V f dx = − ∫ M ⟨ V , ∇ M f ⟩ p dx , (2)where f : M → R is a scalar function and V : M → T M is a vectorfield tangent to the surface. In the following, to simplify the notationwe will drop the subscript M whenever clear from the context. Wefurther denote by L p (M) = { f : M → R | ∫ M | f ( x )| p dx < ∞ , p ≥ } the function space of p -integrable functions on M . The total variation of a differentiable function f on M is defined as TV ( f ) = ∫ M ∥∇ f ( x )∥ dx , (3)where ∥ · ∥ is the norm induced by the inner product ⟨· , ·⟩ p . Thisfunctional quantifies the amount of oscillations in a given signal. Toaccount for non-differentiability of f , we adopt a weaker definitionof TV in terms of continuously differentiable test vector fields U : M → T M [Ben-Artzi and LeFloch 2006]: TV ( f ) = sup U (cid:26)∫ M f ( x ) div U ( x ) dx : ∥ U ( x )∥ ≤ , ∀ x (cid:27) . (4)Alternatively, one can define an anisotropic variant of the standardTV by considering U ( x ) to be bounded in the L ∞ norm. We refer tofunctions f of bounded variation, i.e. TV ( f ) < + ∞ , to be in BV (M) .The TV functional is related to fundamental geometric notionssuch as the curvature of level sets of a function, through the cele-brated co-area formula [Fleming and Rishel 1960]. One particularcase of this relation concerns the total variation of indicator func-tions of closed sets; in this case, the TV corresponds to the perimeterof the set (see [Chambolle et al. 2010, Sec. 1.3] for the Euclideansetting), which generalizes to manifolds in a straight forward way: ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020.
Theorem Let E ⊂ M be a measurable closed subset of M , withsmooth boundary ∂ E , and let χ E be its indicator function. Then: TV ( χ E ) = length ( ∂ E ) . (5) Proof : See Appendix A.1.As we show next, this particular class of indicator functions in-cludes the set of generalized eigenfunctions of the operator associ-ated to the TV functional.
In this work we introduce a new invertible transform for signalsdefined on surfaces, with the following key properties: • The forward transform is non-linear; • The inverse transform is linear; • Any signal f ∈ BV (M) on the surface can be decomposedand reconstructed up to a desired accuracy; • Discontinuities in f are preserved by the transform.The following considerations closely follow the ideas of [Gilboa2013, 2014] in the Euclidean setting, but extend them to functions onmanifolds to enable discontinuity-preserving geometry processing. Consider the following energy for surface signals u , u ∈ BV (M) : E ( u ) = ∫ M ( u ( x ) − u ( x )) dx + tTV ( u ) , t > u ∈ BV (M) E ( u ) . (7)The above seeks to minimize a functional over two terms: a fidelityterm , which quantifies the distance from an initial function u , and a regularization term which measures the total variation of the soughtminimizer. The regularization parameter t balances between thetwo terms. In image processing, where M = S ⊂ R , Eq. (7) isknown as the ROF model for image denoising [Rudin et al. 1992],and dissipates noise in flat regions of an image u while preservingobject contours.Iteratively applying Eq. (7) to the output of a previous solutionwith small increases of t can be seen as a discretization of a timecontinuous total variation flow : ∂ u ∂ t = div (cid:16) ∇ u ∥∇ u ∥ (cid:17) in M , ⟨∇ u , (cid:174) n ⟩ = ∂ M , u ( , x ) = u , (8)where (cid:174) n is the unit vector normal to the boundary ∂ M . Equations (8)describe a physical process of isotropic diffusion within the levelsets of the initial function u , without any diffusion across them (seeFigure 2). In R n , it is guaranteed that the flow converges in finitetime to a constant solution under mild assumptions [Andreu et al.2002, Cor. 1]; no proofs exist for manifolds, but in our experimentswe always observed convergence.To give sense to Eq. (8) for nondifferentiable u and points atwhich ∇ u ( x ) =
0, one needs to turn to the subdifferential , which for t = t = . t = . t = . t = . Fig. 2. Evolution of a disk-shaped function along the TV flow on a sphericalsurface. Since the diffusion only acts inside the level sets, the disk preservesits shape and loses contrast linearly until it vanishes. a convex functional J (TV in our case) over a function space X isformally defined as the set: ∂ J ( u ) = (cid:8) p ∈ X ∗ | J ( v ) − J ( u ) ≥ ⟨ p , v − u ⟩ , ∀ v ∈ X (cid:9) , with X ∗ being the dual space of X . For smooth u with non-zerogradient one obtains: ∂ TV ( u ) = (cid:26) − div (cid:18) ∇ u ∥∇ u ∥ (cid:19)(cid:27) . (9)We refer the reader to [Bellettini et al. 2002] for more details on theTV flow in a Euclidean setting. A generalized eigenfunction of the operator ∂ J is a function u with ∥ u ∥ = λ ∈ R (its corresponding generalized eigenvalue): ∂ J ( u ) ∋ λu , (10)where J = TV in our case. This relation generalizes the classicalcharacterization of eigenpairs of a linear operator: Remark . For the Dirichlet functional D ( u ) = ∫ M ∥∇ u ∥ , oneobtains ∂ D ( u ) = {− ∆ u } , such that the analogy to Eq. (8) becomesthe standard diffusion equation, and Eq. (10) reduces to the standardeigenvalue problem ∆ u = λu (up to a sign). Let us recall the ideas from [Gilboa 2013; Gilboa et al. 2015] to moti-vate a spectral decomposition using the TV flow. For eigenfunctions u : M → R with eigenvalue λ , we get a simple analytic solutionto the flow of Eq. (8): u ( t , x ) = (cid:40) ( − tλ ) u ( x ) for t ≤ λ t ,until it goes to zero. In other words, TV eigenfunctions preservetheir shape (up to rescaling) under the action of the flow; Figure 2illustrates one such eigenfunction.The first derivative of u ( t , x ) w.r.t. time is piecewise-constant: ∂ t u ( t , x ) = (cid:40) − λu ( x ) for t ≤ λ u ( t , x ) vanishes (see Figure 3): ∂ tt u ( t , x ) = (cid:40) + ∞ for t = λ ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020. onlinear Spectral Geometry Processing via the TV Transform • 199:5
Fig. 3. Analytic solution (left) to the 1D flow for a TV eigenfunction u witheigenvalue λ = , first and second derivatives (resp. middle and right). The behavior of the second derivative ∂ tt u ( t ) for an eigenfunction u resembles the notion of spectrum in Fourier analysis, since sinu-soidal functions (eigenfunctions of the Laplace operator) correspondto deltas in the frequency domain. This motivates the definition: Definition We refer to the distribution over time ϕ t = ∂ tt u ( t ) as the spectral representation of u , and to each individual distributionat the time scale t as a spectral component . As we have seen, for a TV eigenfunction with eigenvalue λ , thespectral representation has a single component: a delta at t = λ .However, not every spectral component ϕ t is guaranteed to identifya TV eigenfunction. Further, at any given time scale t , the spectralcomponent ϕ t : M → R still has a spatial extent over M . Toshow the entire spectral activity across different time scales, we canintegrate each ϕ t and simply plot the function: s ( t ) = ∫ M | ϕ t ( x )| dx . (11)Looking at the behavior of the function s ( t ) can be rather useful inpractice, since it helps in defining the filtering operations requiredby the applications; examples of s ( t ) are shown in Figure 1.The spectral representation can be translated back to the primaldomain by the linear operation of integrating over time (i.e., acrossall spectral components ϕ t ) and adding the mean of u : u = ∫ ∞ tϕ t dt + mean ( u ) . (12)We refer to [Gilboa 2013, 2014] for more details (with the slightdifference that we decided to define the spectral decompositionwithout the factor t and include it in the reconstruction instead). In light of the previous discussion that TV eigenfunctions are thefundamental atoms of the spectral representation of Def. 1 (in thesense that they yield peaks in the spectral domain), one might expectto represent arbitrary signals as linear combinations of such atoms.Although there seems to be a vast number of TV eigenfunctionsto allow such a representation , computing a sparse representationturns out to be difficult since, e.g., a Rayleigh principle to computeorthogonal eigenfunctions does not hold [Benning and Burger 2013].Interestingly, Def. 1 allows to show that the gradient flow withrespect to a regularization with suitable properties does yield sucha decomposition with the ϕ t representing differences of eigenfunc-tions, see [Burger et al. 2016]. While the TV in Euclidean spaces See e.g. [Steidl et al. 2004] for a proof that the Haar basis is a proper subset of theTV eigenfunctions on the real line, or [Alter et al. 2005] for a characterization of setswhose characteristic function is a TV eigenfunction.
Fig. 4. Spectral TV transform on a rubber ducky. The original signal (topleft) is composed of three indicator functions of different regions, repre-sented in the TV spectral domain as three separated peaks (top right). Eachplot visualizes the spectral activity of the corresponding surface signal. Byfiltering the peaks in the spectral domain we can single out each componentat a different scale, obtaining one region at a time (bottom row). of dimension larger than 1 does not satisfy the required properties,the spectral representation ϕ t has been shown to still yield a mean-ingful, edge-preserving, data-dependent multi-scale representationin the field of image processing [Gilboa 2014; Gilboa et al. 2015],which is why we utilize it on manifolds in a similar fashion. To compute the spectral representation ϕ t of a given function f , weuse an implicit Euler discretization of the TV flow: We compute anew iterate by solving Eq. (7) using the previous iterate in the datafidelity term starting from u ≡ f . The procedure is summarized inAlgorithm 1. Algorithm 1:
Spectral TV decomposition input : signal u , step length α , no. of time steps N output : spectral representation ϕ t begin for t = N do u ( t ) Alд . ←−−−−− min u TV ( u ) + α || u − u || output ϕ t ←− ∂ tt u ( t ) increase α update u ←− u ( t ) end end The algorithm simply evolves the input signal by N discretesteps along the TV flow; each iteration moves a step forward, withdiffusion time equal to α . As changes happen quickly for small t andtend to become slower for larger t we iteratively increase the stepsize α of the evolution. Subsequently, the spectral representation ϕ t is constructed incrementally using finite differences, such that theintegral of Eq. (12) becomes a simple weighted sum over the ϕ t . Wegive selection strategies for α and N in Appendix B.4. Remark . This yields a multiscale representation for the input,since each diffusion time t captures a feature at a different scale – ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020. namely, the feature that vanishes at t . Here, time scales play the samerole as a wavelength in classical Fourier analysis. See Figures 4 and 5 for examples of multiscale TV decompositionof color signals on two different surfaces.
Algorithmic considerations . Algorithm 1 is unconditionally sta-ble, and line 3 corresponds to taking an implicit Euler step along thediscretized TV flow. This property allows us to use arbitrary valuesfor the step size α (therefore, the resolution of the decomposition)without compromising convergence. In practice, in our experimentswe use a variant of Algorithm 1 and study the evolution of the flowbackwards in time by following a nonlinear inverse scale approach,that starts from a constant solution and converges to the initialfunction u . The spectral decompositions of both approaches arevery similar in practice and have been proven to be identical un-der certain assumptions, see [Burger et al. 2016]. The inverse scalevariant comes with additional numerical robustness. The completealgorithm, along with a detailed discussion on the two flows and onthe role of α are given in Appendices B.1 and B.2. We conclude by giving a more complete characterization of theTV eigenfunctions on surfaces. On the image domain, Bellettini etal. [2002] proved that any indicator function χ S : R → { , } of a convex closed subset S ⊂ R is a TV eigenfunction with eigenvalue λ = Per ( S )| S | if it satisfies the isoperimetric inequality:max q ∈ ∂ S k ( q ) ≤ Per ( S )| S | , requiring the maximal curvature k of ∂ S to be less than or equal tothe perimeter-area ratio of S , intuitively excluding oblong regions.Such results have later been generalized to R n as well as tounions of such convex sets, where the convex sets need to besufficiently far apart for the TV spectral decomposition to yielda single δ -peak for each of the subsets, see [Alter et al. 2005].We conjecture (and consistently con-firmed in all our experiments) thatthis result can be generalized to sur-faces by replacing convex subsets with geodesically convex subsets ofsurface M , making the examples in the inset figure eigenfunctionsof the TV on a manifold. A patch C ⊂ M is geodesically convexif, given any two points p , q ∈ C , there exists a unique minimum-length geodesic connecting them which lies entirely in C . Therefore source filtered A filtered B Fig. 5. The source surface signal (left) is transformed into a spectral TVrepresentation, and processed with two band-pass filters to obtain a cowwithout eyes and nostrils (A) and without spots (B).
Input
Fig. 6. A U-shaped function (left) is decomposed as a linear combinationof TV eigenfunctions corresponding to two geodesically convex sets, mani-fested as two peaks in the TV spectrum. Functions are color-coded, growingfrom blue (negative values) to red (positive values). in certain practical cases (e.g., the cow’s nostrils in Figure 5) oneactually gets TV eigenfunctions in the spectral decomposition.
Input ϕ ϕ ϕ ϕ Fig. 7. Spectral TV decomposition of the mean curvature function on thebimba (left). The spectral components ϕ t are visualized at increasing index.The curve plot shows the spectral activity function s ( t ) of Eq. (11) . Low fre-quency components localize on semantically relevant regions (the chignon),while higher frequencies localize along feature lines (the hair strands). In Figure 6 we show an example of indicator function of a non -geodesically convex set (U-shaped, left) that can be represented by alinear combination of indicator functions of geodesically convex sets(right). In Figure 7 we show the spectral components of a more com-plex signal (the mean curvature); in this case, the components arenot geodesically convex but we still get a meaningful decompositionthat can be used for mesh processing.
We implemented the framework in Matlab 2019a, and ran all the ex-periments on an Intel i7-4558U CPU with 16 GB RAM. A discussionon runtime performance is provided in Appendix C. v j v k v i e ik e ij We discretize shapes as manifold triangle meshes ( V , E , F ) , sampled at vertices v i , i = { ... | V |} , withthe constraint that each pair of triangles ( t i , t j ) ∈ F shares at most one edge e ∈ E . Functions f : M → R are approximated at vertices with piecewise-linear basis elements; ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020. onlinear Spectral Geometry Processing via the TV Transform • 199:7 tangent vector fields U : M → T M are constant within eachtriangle.The gradient operator G ∈ R | F |×| V | takes a function value pervertex v i and returns a 3D vector per triangle t ijk ∈ F : G f ( t ijk ) = (cid:18) v ⊤ j − v ⊤ i v ⊤ k − v ⊤ i (cid:19) ⊤ (cid:18) ∥ e ij ∥ ⟨ e ij , e ik ⟩⟨ e ij , e ik ⟩ ∥ e ik ∥ (cid:19) (cid:18) f ( v j ) − f ( v i ) f ( v k ) − f ( v i ) (cid:19) where the notation is as per the inset figure. The divergence operator D ∈ R | V |× | F | takes a 3D vector per triangle and returns a valueper vertex. As the adjoint of the gradient, it is discretized as : D = − A − G ⊤ T , (13)where A is a diagonal matrix of local area elements at each vertex(shaded area in the inset) and T is a diagonal matrix of triangle areas.These discretizations follow the standard formulations as described,e.g., in [Botsch et al. 2010].To define these operators on point clouds, we base our approachon the construction of a graph (V , E) , built locally at each pointin a k -neighborhood. The edges of the graph are weighted by aGaussian kernel w ij = e − ∥ xi − xj ∥ δ , for an appropriate δ dependingon the point set density. Functions f : V → R (E) are defined onthe nodes of the graph while vector fields F : E → R (V) are definedon the edges. The gradient and divergence operators are defined as: G ij = w ij ( f i − f j ) ∀ ( i , j ) ∈ E (14) D ( F i ) = (cid:213) j ∈( i , j )∈E w ij ( F ij − F ji ) (15)In order to deal with low-quality meshes, we adopted the ap-proach of [Sharp et al. 2019]; a more detailed discussion with anexperimental comparison is given in Appendix B.4. The energy to minimize in line 3 of Alg. 1 is convex, but non-differentiable. To efficiently solve it, we discretize the problem anduse the primal-dual hybrid gradient algorithm [Chambolle and Pock2011], which is well-suited for problems of this form. The optimiza-tion procedure for an input signal u is summarized in Algorithm 2:the proximal operator prox F ∗ is the projection of each component Algorithm 2:
Primal-Dual Hybrid Gradient (PDHG) input : u , q , σ , τ > θ ∈ [ , ] output : minimizer u ∗ of Eq. (7) begin while not converged do u k + = prox G , τ ( u k − τ Dq k ) ¯ u k + = u k + + θ ( u k + − u k ) ; q k + = prox F ∗ ( q k + σ G ¯ u k + ) end return u ∗ = u k + end of its input onto the unit L ball: prox F ∗ ( q j , : ) = (cid:40) q j , : ∥ q j , : ∥ ≤ q j , : ∥ q j , : ∥ otherwise , (16) Fig. 8. Diffusing vertex coordinates on the cat model. An anisotropic TVpenalty will promote piecewise-flat solutions, where oscillations are min-imized along the coordinate axes. Therefore, one gets different results byrotating the global reference frame; for the source model on the left, we showsolutions under two different axis rotations (middle and right respectively). where each q j , : ∈ R is a vector at triangle j . The proximal operator prox G , τ has a closed-form solution given by: prox G , τ ( u ) = u + τα u + τα . (17)We initialize the algorithm with u coming from the previousiteration in Algorithm 1 and q = , and stop the iterations whenthe absolute change in energy is below a small threshold. We set θ = .
5, which appeared optimal in our numerical experiments.
We now study the evolution of the surface M along the TV flow. Adirect way to do this is to evolve the ( x , y , z ) vertex coordinates of M ; the same approach is also used to define, e.g., the mean curvatureflow (MCF) and its variants [Desbrun et al. 1999; Kazhdan et al. 2012;Taubin 1995]. This approach, however, gives rise to two main issues.First, diffusing the vertex coordinates naturally affects the metricof M at each diffusion step; this, in turn, modifies the functional TV ( f ) = ∫ M ∥∇ M f ∥ , since the integration domain M and theattached intrinsic operators now vary along the flow. Recomputingthe metric-dependent operators at each iteration of the algorithm isa possible solution, but can be highly inefficient.Secondly, one must take care of the choice of the TV regularizationfor ( x , y , z ) . In fact, for f being the coordinate function on the man-ifold, ∇ f ( x ) becomes a tensor for which a suitable norm ∥∇ f ( x )∥ in the integral of the total variation needs to be defined. Dependingon such a choice one can encourage different types of collaborativegradient sparsities [Duran et al. 2016]. For example, by choosing afully anisotropic TV (i.e. summing the absolute value of all entriesin ∇ f ( x ) ), the transform will be orientation-dependent: the diffusedvertex coordinates will tend to align to the global reference frameof R , giving different results depending on the orientation of theinitial shape; see Figure 8 for examples. While this latter point canbe exploited for certain stylization tasks (as we show in Section 7),it might be an undesirable effect in many other applications. In this paper we advocate the use of normal fields for encoding thegeometry. This choice has three main benefits: (1) the domain uponwhich the operators are defined remains fixed throughout the flow;(2) using normals endows us with rotation invariance ; and (3) we getinteresting connections with other existing flows. Normal fields inconnection with the TV functional have been considered before, e.g.in [Zhang et al. 2015; Zhong et al. 2018]. However, these approaches
ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020.
Fig. 9. Our flow induces the formation of clusters in the normal field. Thetop row shows the original surface, the bottom row is the filtered result.
Middle : the local density of the Gauss map (heatmap growing from white tored) illustrates the formation of clusters in normal space.
Right : the coverageof S (front and back sides visualized) shows that the diffused normals areless spread across the sphere than in the original surface. do not consider or analyze the geometric flow as we do here, and donot provide a natural and interpretable multiscale decomposition. Γ −→ Our input signal is the normal vectorfield (cid:174) n over M . The field (cid:174) n is iden-tified by the image of the Gauss map Γ : M → S , mapping points on thesurface to points on the unit 2-sphere(see inset). The corresponding time continuous flow is: (cid:40) ∂ Γ ∂ t = div S (cid:16) d Γ ∥ d Γ ∥ (cid:17) in S , Γ ( , x ) = (cid:174) n , (18)where d Γ is the differential of the Gauss map and div S is definedwith the metric of M , but applied to tangent fields on S . By dif-fusing normals according to Eq. (18), the point distribution over S changes to form well-separated clusters. Intuitively, the flow tendsto align the directions of normal vectors to favor piecewise-constantsolutions; see Figure 9 for an example. To obtain the filtered em-bedding for M , the diffused normal field is integrated as describedfurther in this section.Interpreting this approach from a differential geometric perspec-tive provides us with a useful characterization of the normal TVflow in terms of surface curvature. In particular, the TV functionalfor the normal field (cid:174) n reads: TV ((cid:174) n ) = ∫ M ∥ d Γ ( x )∥ dx = ∫ M (cid:113) k + k dx , (19)where k , k are the principal curvatures at each point. The lastequality follows from the fact that k , k are the eigenvalues of d Γ .As we show below, the final expression in Eq. (19) resembles thefunctional related to other known flows in the literature. We alsosee from the same expression that the TV energy of the normal fielddoes not depend on the orientation of the normals in R ; therefore,the resulting flow is invariant to global rotations of M . source Fig. 10. Comparison between MCF (top row) and normal TV flow (bottomrow) on a 1D closed contour of a dog. The circles visualize the image of theGauss map on S for each shape; the color code is defined as on the sourceshape (leftmost column). The MCF equally spreads the normals, eventuallyconverging to a circle, while the TV flow forms well-separated clusters (the5 points on the circle), converging to a piecewise-flat shape. Implementation.
By working with normal fields (cid:174) n we are movingto manifold- valued functions (the manifold being the unit sphere).In our tests we treat each (cid:174) n ( x ) as an element of R , and add anormalization step projecting (cid:174) n ( x ) back onto S at each iteration.In the discrete setting, the normal vector field n ∈ R | F |× isconstant within each triangle. Its gradient (needed in Algorithm 2,line 5) is understood as the channel-wise jump discontinuity overeach edge in each triangle, and the associated discrete operator canbe assembled as a | E | × | F | matrix; see Appendix B.3 for details. Recovering vertex positions.
To recover a new embedding fromthe filtered normals n , we follow [Prada and Kazhdan 2015]:min v ∈ R | V |× ∥ Gv − w ∥ F + ϵ ∥ v − v ∥ F , (20)where v are the old vertex coordinates, and w = Gv −⟨ Gv , n ⟩ r ⊙ n is the orthogonal component of Gv with respect to n ( ⟨· , ·⟩ r and ⊙ operate row-wise). Eq. (20) seeks the set of vertices whose gradientmatches the vector field w , while at the same time staying closeto the initial vertices. In our tests we set ϵ = − to allow largedeviations from the initial shape. Solving Eq. (20) boils down tosolving a screened Poisson equation, giving the sparse linear system: v = ( ϵ A − L ) − ( ϵ Av − Dw ) . (21)Here, L = − DG corresponds to the standard linear FEM discretiza-tion of the Laplace-Beltrami operator. We compare our flow to other geometric flows used in graphics.
Mean curvature flow.
The normal TV flow is closely related tothe MCF, which describes the evolution of a surface under theminimization of the area functional (membrane energy): E A ( u ) = ∫ M dx , (22)where u is now the ( x , y , z ) embedding in R of the surface M , andwe assume M has no boundary. The relation to our flow emergesby approximating the membrane energy with the Dirichlet energy ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020. onlinear Spectral Geometry Processing via the TV Transform • 199:9
Fig. 11. Normal TV flow on the bimba. Sharp geometric features getsmoothed in a “vertical” fashion, preserving the overall structure. For exam-ple, the chin and pony tail maintain their shape throughout, while the braidand nose slowly disappear without affecting the neighboring geometry. E A ( u ) ∝ D ( u ) = ∫ M ∥∇ u ∥ , whose gradient flow is the PDE: ∂ u ∂ t = ∆ u . (23)Compared to Eq. (8), this PDE uses the linear operator ∆ = div ∇ instead of the nonlinear operator div ∇∥∇∥ .Denoting with H = k + k the mean curvature, and observingthat ∆ u = H (cid:174) n , the above PDE corresponds to surface motion in thenormal direction, with speed proportional to curvature H . Criticalpoints of the functional are minimal surfaces, which for genus-0convex shapes correspond to spheres.In comparison, surface evolution along the normal TV flow col-lapses the surface to a flat shape, with the direction of motion per-pendicular to the level sets of the normal vector field (cid:174) n . By lookingat the effect of the flow on (cid:174) n , the key difference is that the MCFtends to move points on S so as to keep them equally spaced,while the normal TV flow generates well-separated point clusters.We illustrate the comparison in 2D in Figure 10; an example of ourfeature-preserving flow in 3D is given in Figure 11. Willmore flow.
Perhaps more interestingly, we observe that ourfunctional TV ((cid:174) n ) resembles the Willmore (or thin plate) energy[Bobenko and Schröder 2005], a higher-order version of the areafunctional, defined as: E W ( u ) = ∫ M H dx = ∫ M ( k + k ) dx − (cid:24)(cid:24)(cid:24)(cid:24) π χ (M) , (24)where again u encodes the coordinate functions in R . The lastterm is the constant Euler characteristic of M , and thus it doesnot affect the flow. Compared to Eq. (19), we see that the integrandonly changes by a square root operation. In this sense, we interpretthe normal TV flow as a sparsity-promoting variant (in the normaldomain) of the Willmore flow.The linearization of Eq. (24) leads to a fourth-order bi-Laplaciandiffusion equation: ∂ u ∂ t = ∆ u . (25)Our problem also involves a fourth-order PDE with respect to thecoordinates; however, it is hard to solve directly due to the nonlin-earity of our operator. By phrasing the PDE in terms of the normalsinstead of the coordinates, we decouple the problem: first we inte-grate the TV flow in the space of normals, which corresponds to the Flow Functional PDE
MCF E A ( u ) = ∫ dx ∂ u ∂ t = ∆ u Willmore E W ( u ) = ∫ ( k + k ) ∂ u ∂ t = ∆ u TV E TV ((cid:174) n ) = ∫ (cid:113) k + k dx ∂ Γ ∂ t = div S (cid:16) d Γ ∥ d Γ ∥ (cid:17) Table 1. Comparison between MCF, Willmore, and our normal TV flow interms of functional involved and the associated PDE. T V W i ll m o r e M C F convergence Fig. 12. Evolution of the bunny surface along different geometric flows.Top to bottom: conformalized mean curvature flow [Kazhdan et al. 2012],conformal Willmore flow [Crane et al. 2013], and our normal TV flow. TheTV flow differs from the others as it preserves the underlying structure, anddetails are removed gradually. At convergence, we get a flat surface sincethe image of the Gauss map is clustered into a single point on S . second -order PDE of Eq. (18), and then we integrate again to recoverthe vertex coordinates via the screened Poisson equation.In Table 1 we summarize the three flows. A qualitative comparisonis shown in Figure 12, where we use stable variants of the existingflows that prevent singularities. We remark here that both the Laplace-Beltrami operator ∆ M andthe subdifferential ∂ TV explored in this paper are special cases of amore general parametric operator, the p -Laplacian [Lindqvist 2006]: ∆ p = div (∥∇∥ ( p − ) ∇) , (26)with the associated p -Dirichlet energy: E p ( u ) = p ∫ M ∥∇ u ∥ p dx . (27)In particular, the p -Dirichlet energy with p = p = p -Laplacian flow of a unit ball, we observe convergence to a cubefor p = p = p ≫
2, which is a consistent behavior withthe definition of p -norms (see Figure 13). From the example we seeclearly that the TV flow tends to piecewise-flat solutions, and furtherobserve that other p -Laplacian flows carry other properties that maybe worth investigating in the future. We refer to [Bühler and Hein ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020. p = p = p = p = . MCF TV
Fig. 13. Evolution of the vertex coordinates of the unit L sphere along p -Laplacian flows, for different values of p . The flows are simulated withthe explicit Euler method. For p = the sphere remains a sphere, while itbecomes a piecewise-flat surface with the TV flow ( p = ). Method TV Energy
Source 43 . Ours . [Sun et al. 2007] 27 . . . . . Table 2. Quantitative comparison between our TV spectral filtering methodand the feature-preserving approaches from Figure 14. For each method wereport the TV energy of Eq. (19) , which quantifies the amount of oscillationsin the geometry after filtering; larger values correspond to spurious noise. p -Laplacian and to[Cohen and Gilboa 2019] for an approach defining a p -Laplacian-based spectral decomposition in Euclidean spaces. The TV spectral framework lends itself well to a variety of applica-tions. In these, we follow the same computational procedure:(1) Compute the forward TV transform of a given signal (possibly,the geometry itself);(2) Apply a filter to the resulting spectral components;(3) Compute the inverse TV transform.Letting u : M → R k be a k -channel signal on surface M , ϕ t beits spectral representation, and I be a (possibly nonlinear) filteringoperator, the filtered signal ˜ u is computed as:˜ u = ∫ tI ( ϕ t ) dt . (28)Note that the above formula is meant to be generic, while in practicethe type of TV regularization used to construct ϕ t plays an importantrole as we shall see below.We showcase our method on a range of relevant applications forgraphics; for each application, we compare with other approachesfrom the state of the art to better position our method in terms ofachievable quality. While these results are intended to demonstratethe flexibility of the framework, they are by no means exhaustiveand could serve as a basis for follow-up work. Perhaps the most straightforward application is the direct filteringof signals and geometric features living on the given surface. Onesuch example is shown in Figure 1, where we apply a band-pass “square window” filter of the form: I ba ( ϕ t ) = (cid:40) ϕ t t ∈ [ a , b ] [ a , b ] is the range of time scales left untouched by the filter;the normal field (cid:174) n is decomposed by applying Algorithm 2 to eachcomponent of (cid:174) n separately and coupling the channels by projectingon the unit sphere at each step. The above filter relies upon one mainproperty of the TV decomposition: The separation of geometricdetails into different spectral bands, ensured by the underlyingassumption that most interesting shapes have a sparse gradient. Thissimple filter already provides high-quality selective suppressionof semantic detail while preserving sharp geometric features. InFigure 14 we compare to five other top performing methods fromthe literature,specifically tailored to edge-preserving filtering, whilewe show a comparison with Laplacian-based filtering in Figure 15. A remarkable application of the spectral TV framework is to transfersharp geometric details between shapes. Consider two surfaces M , N and a bijection π : M → N between them. Further, let ϕ t and ψ t be the TV spectral representations of the normal fields on M and N respectively. The idea is to transfer some componentsof ϕ t to ψ t via π , and then compute the inverse transform of ψ t toobtain a new version of shape N with additional semantic details.For this to work, recall that ϕ t : M → R is a vector field on M for any given t ∈ R + , and similarly for ψ t : N → R . Therefore, wecan transfer details from M to N by evaluating the integral:˜ u ( x ) = ∫ t ( ψ t ( x ) + I ba ( ϕ t ) ◦ π − ( x )) dt (30) = u ( x ) + ∫ ba t ( ϕ t ◦ π − ( x )) dt , (31)where u is the original normal field on N , ˜ u is the newly synthesizednormal field, and I ba is a band-pass filter that selects the desiredfeatures to borrow from M . The integral in Eq. (30) is a directapplication of the reconstruction formula of Eq. (12) to an additivelyupdated spectral representation ψ t . The key feature of this procedureis that, despite its simplicity, it enables a form of selective transferof details. Importantly, the details must not be necessarily localizedon the surface, as long as they are well localized in scale space.In Figure 16 we show an example of detail transfer. In Figure 17we compare our approach to the Laplacian-based detail transferapproach of [Sorkine et al. 2004], where differences between Lapla-cian coordinates (i.e. differences between normals scaled by meancurvature) are transferred between shapes. For these tests, the map π was computed using the method of [Melzi et al. 2019] startingfrom 15 hand-picked landmark correspondences. We further exploit some basic properties of the TV flow for artisticrendition of 3D shapes. In particular, we show how with a properchoice of the TV functional , we get the desired effect of mimickingthe stylistic features typical of voxel art. A similar effect of “cubicstylization” was achieved in [Huang et al. 2014] for constructing
ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020. onlinear Spectral Geometry Processing via the TV Transform • 199:11
Source
Ours [Zheng et al. 2011] [Fleishman et al. 2003] [Sun et al. 2007] [Zhang et al. 2015] [Jones et al. 2003]
Fig. 14.
Top row:
Example of high-frequency detail removal, implemented as a low-pass filter on the spectral TV representation of the source shape (leftmost).We get a “rejuvenation” effect on the face model (second column), and compare our solution with state-of-the-art approaches, specifically tailored foredge-preserving mesh filtering (third to fifth column).
Bottom row:
To better highlight the differences, we plot the mean curvature (clamped to [− , ] andcolor-coded, growing from blue to red) on the filtered surfaces. Our method is the most effective at removing wrinkles and obtains sharper feature linescompared to other approaches, which either tend to over-smooth or retain spurious vertex noise. A quantitative comparison is reported in Table 2. Laplacian TV
Fig. 15. Comparison between low-pass filtering of a source shape (left)using 300 Laplacian eigenfunctions (middle), and with our spectral TVrepresentation using just
11 spectral components out of 25 (right). The fullshape has approximately k vertices. Laplacian-based smoothing correctlyremoves details, but also has a smoothing effect on the rest of the shape. Incontrast, our approach correctly preserves the underlying structure. Source Target Result
Fig. 16. Example of detail transfer. The wrinkles on the source shape X (left)are captured with a high-pass filter in the spectral TV decomposition of X , added to the spectral representation of the target (middle), and finallyresynthesized (right). The TV decomposition allows to inject tiny geometric details into the target, without affecting the other bands (e.g., nose, mouthand ears are not affected by the transfer). Note that source and target shapedo not have the same mesh connectivity, nor the same number of vertices. polycube maps for texturing, and in [Liu and Jacobson 2019] formodelling purposes. These approaches minimize an as-rigid-as-possible deformation energy with L regularization on the normalfield. Similarly to these methods, our modified meshes retain anyattribute defined on the original shape (e.g. UV coordinates fortextures). Further, since in our case the stylization is the result ofan evolution process, this can be stopped at any time to attain thedesired sharpness level. A simple way to get cubic stylization isto consider the ( x , y , z ) coordinates as separate scalar functionsregularized via the sum of anisotropic TV penalties. This causes thesurface to flatten during the diffusion, and eventually shrink to asingle point at convergence. Due to the strict dependence on theglobal reference frame, different orientations of the surface in R will give different results, as in [Huang et al. 2014; Liu and Jacobson2019]. Shrinkage can be avoided by diffusing the normal vector fieldinstead of the vertex positions. This way, the flow only modifies thenormal directions, and solving the screened Poisson equation (21)approximately preserves the volume. See Figures 17 and 19 for someof our stylization examples. Finally, we demonstrate the application of the spectral TV frame-work to more classical tasks of geometry processing, namely surfacedenoising and enhancement. For these tests, we use both trianglemeshes and point clouds. Once the spectral decomposition is com-puted, filtering and reconstruction can be performed in real time,as shown in the accompanying video.
Fairing.
When dealing with point clouds, denoising and fairingtasks are notoriously challenging due to the difficulty of distinguish-ing feature points from outlier noise. This confusion often requirespre-processing steps to localize sharp features and filter out therest. In our setting, the edge-aware property of the TV functional
ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020.
Source[Sorkine et al. 2004]
Ours
Target
Fig. 17.
Replacing the inscription from a source (top left) to a target model(top right), using Laplacian coating transfer [Sorkine et al. 2004] (bottomleft), and our detail transfer approach (bottom right). The smoothing appliedto the target causes the former approach to fail at preserving edges, whilethe transfer of high frequencies in the TV spectrum preserves them. Thetarget model was obtained by applying our cubic stylization method ofSection 7.3 to the source model before superimposing the inscriptions. allows us to automatically damp oscillations due to noise, while atthe same time retaining the true details of the underlying surface.We formulate denoising as a simple low-pass filter in the spectralTV domain. The input signal is the normal field (cid:174) n : M → R ; when M is discretized as a point cloud, normals can be estimated via totalleast squares [Mitra and Nguyen 2003]. See Figure 18 for an exampleof denoising of point clouds, where we compare with the recentbilateral filtering approach of [Digne and de Franchis 2017]. Enhancement.
While denoising can be seen as a suppression ofhigh-frequency components in the spectral representation of (cid:174) n ,one can similarly obtain feature enhancement by amplifying thosecomponents; see Figure 20 for an example. Here we compare to thegradient-based mesh processing approach of [Chuang et al. 2016],which implements feature sharpening and smoothing by scaling thegradient of coordinate functions and solving for new coordinates viaa Poisson equation. Conceptually, our spectral TV representationallows more accurate control on specific features, since these arewell separated in the spectrum; differently, the Laplacian-basedapproach adopted in [Chuang et al. 2016] encodes geometric featuresby spreading them out across the entire spectrum. Localized operations.
We conclude by mentioning that our frame-work also supports localized operations by means of masks definedover regions of interest, which could be useful in view of an inter-active application for mesh editing (see the supplementary video).Let R ⊂ M be a region of M , and let I be an arbitrary spectral filter.Then, for a signal u with spectral components ϕ t , localized spectralfiltering can be carried out simply by evaluating:˜ u ( x ) = ρ ( x ) ∫ M tI ( ϕ t ( x )) dt + ( − ρ ( x )) u ( x ) , (32)where ρ : M → [ , ] is an indicator function for R . input 0.059 Ours 0.022 [Digne and de Franchis 2017] 0.027ground truth
Fig. 18. The oriented point cloud on the top left (obtained by additiveGaussian noise on the vertex coordinates) is denoised by applying a low-pass filter to the spectral TV representation of its normal field (bottom left).On the bottom right we compare to Digne [2017]. At the top-right of eachpoint cloud, we plot its normal field on the unit sphere. Points are color-coded according to their cosine similarity (light to dark red, with valuesreported in the figure) with the ground truth normals (in solid cyan)
We presented a new spectral processing framework for surface sig-nals and 3D geometry. The framework is based on ideas developedin the last decade within the nonlinear image processing commu-nity, but whose extension and application to geometric data havebeen lagging behind to date. The main contribution of this paper isto show how these ideas can serve as valuable tools for computergraphics, and to relate the resulting observations and algorithmsto more well known paradigms such as Laplacian-based spectralapproaches and shape fairing via geometric flows. Probably themost exciting aspect of the framework lies in its interpretability interms of spectral components with an intuitive canonical ordering.Together with the jump-preserving property, one gets the forma-tion of frequency bands that tend to cluster together features withsimilar semantics. For example, the hair strands on a head modelwill be captured in a higher-frequency band than the ears; both willbe captured sharply. The general approach to mesh filtering that wedescribed also retraces the steps of classical pipelines adopted in thelinear case (transform – filter – resynthesize), making our nonlinearframework more accessible and of potentially broader impact.
From a practical perspective, perhaps our weakest point is the needto define a schedule for the parameter α , which determines theresolution of the spectral representation, influencing the design offilters. While in most cases an exact tuning is not required, it mightbe possible to define an adaptive scheme to automatically tune thescale of α to the feature landscape of the signal. Similarly, a moresophisticated filter design to achieve complex effects is possible with ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020. onlinear Spectral Geometry Processing via the TV Transform • 199:13
Fig. 19. Example of cubic stylization. The scene on the left is evolved along the anisotropic normal TV flow (i.e., each dimension of the normal vectors istreated as a separate scalar function, and these are regularized jointly via the sum of anisotropic TV penalties), which forces the normals to align with theglobal coordinate axes; the final result is shown on the right. Each shape was processed individually, and then placed in the scene.
Source
Ours [Chuang et al. 2016] [Chuang et al. 2016] w/ smoothing
Fig. 20.
Top row:
Spectral TV enhancement of a lion-shaped vase. Assigninga larger weight to the high-frequency spectral components has an effectof pronouncing the details of the 3D model, while preserving the sharpfeatures.
Bottom row:
Comparison with [Chuang et al. 2016]. Their methodcomes in two variants. On the left, their direct sharpening results in a noisymesh. On the right, feature enhancement is alternated with smoothing 3times, at the cost of detail loss (e.g. around the whiskers). the adoption of learning-based techniques to predict task-specificfilters , as done in [Moeller et al. 2015] for image denoising. Moregenerally, the integration within geometric deep learning pipelinesis a promising direction that we are eager to pursue. We also an-ticipate that the study of the more general p -Laplacian operatoron surfaces might lead to a richer class of algorithms for geometryprocessing. One promising idea is to define the operator adaptively,selecting a value for p depending, e.g., on curvature. Finally, a morein-depth analysis of the normal TV flow as a robust alternative formesh fairing is essential, due to its feature preservation propertiesand its connection to other geometric flows. ACKNOWLEDGMENTS
The authors gratefully acknowledge the anonymous reviewers forthe thoughtful remarks, and Matteo Ottoni, Giacomo Nazzaro andFabio Pellacini for the artistic and technical support. MF and ER aresupported by the ERC Starting Grant No. 802554 (SPECGEO) andthe MIUR under grant “Dipartimenti di eccellenza 2018-2022” of theDepartment of Computer Science of Sapienza University.
REFERENCES
Yonathan Aflalo, Haim Brezis, and Ron Kimmel. 2014. On the Optimality of Shape andData Representation in the Spectral Domain.
SIAM J. Imaging Sciences R n . Math. Ann.
332 (2005), 239–366.Fuensanta Andreu, Vicente Caselles, Jesú Díaz, and José Mazón. 2002. Some QualitativeProperties for the Total Variation Flow.
J. Funct. Anal.
Computer Vision -ECCV 2014 Workshops . Springer International Publishing, 299–312.Chandrajit L. Bajaj and Guoliang Xu. 2003. Anisotropic Diffusion of Surfaces andFunctions on Surfaces.
ACM TOG
22, 1 (Jan. 2003), 4âĂŞ32.Giovanni Bellettini, Vicente Caselles, and Matteo Novaga. 2002. The Total VariationFlow in RN.
Journal of Differential Equations
Annales de l’Institut HenriPoincare (C) Non Linear Analysis
24 (12 2006), 989–1008.Martin Benning and Martin Burger. 2013. Ground states and singular vectors of convexvariational regularization methods.
MAA
20, 4 (2013), 295âĂŞ334.Alexander I. Bobenko and Peter Schröder. 2005. Discrete Willmore Flow. In
Proc. SGP .Eurographics Association, 101âĂŞes.Nicolas Bonneel, David Coeurjolly, Pierre Gueth, and Jacques-Olivier Lachaud. 2018.Mumford-Shah Mesh Processing using the Ambrosio-Tortorelli Functional.
ComputGraph Forum
37, 10 (Oct. 2018).Mario Botsch, Leif Kobbelt, Mark Pauly, Pierre Alliez, and Bruno Lévy. 2010.
Polygonmesh processing . CRC press.Christopher Brandt, Leonardo Scandolo, Elmar Eisemann, and Klaus Hildebrandt. 2017.Spectral processing of tangential vector fields.
Comput Graph Forum
36, 6 (2017),338–353.Thomas Bühler and Matthias Hein. 2009. Spectral Clustering Based on the GraphP-Laplacian. In
Proc. ICML . 81âĂŞ88.Martin Burger, Guy Gilboa, Michael Moeller, Lina Eckardt, and Daniel Cremers. 2016.Spectral decompositions using one-homogeneous functionals.
SIAM Journal onImaging Sciences
9, 3 (2016), 1374–1408.Antonin Chambolle, Vicent Caselles, Daniel Cremers, Matteo Novaga, and ThomasPock. 2010. An introduction to total variation for image analysis.
Theoreticalfoundations and numerical methods for sparse recovery
9, 263-340 (2010), 227.Antonin Chambolle and Thomas Pock. 2011. A first-order primal-dual algorithm forconvex problems with applications to imaging.
JMIV
40, 1 (2011), 120–145.Ming Chuang, Szymon Rusinkiewicz, and Michael Kazhdan. 2016. Gradient-DomainProcessing of Meshes.
J. Comput. Grap. Tech.
5, 4 (Dec. 2016), 44–55.ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020.
Ulrich Clarenz, Udo Diewald, and Martin Rumpf. 2000. Anisotropic Geometric Diffusionin Surface Processing. In
Proc. VIS . IEEE, 397âĂŞ405.Ido Cohen and Guy Gilboa. 2019. Introducing the p-Laplacian Spectra. (2019).Dorin Comaniciu and Peter Meer. 2002. Mean Shift: A Robust Approach Toward FeatureSpace Analysis.
IEEE PAMI
24, 5 (May 2002), 603–619.Keenan Crane, Ulrich Pinkall, and Peter Schröder. 2013. Robust Fairing via ConformalCurvature Flow.
ACM TOG
32 (2013).Mathieu Desbrun, Mark Meyer, Peter Schröder, and Alan H. Barr. 1999. Implicit Fairingof Irregular Meshes Using Diffusion and Curvature Flow. In
Proc. SIGGRAPH . ACMPress, 317–324.Mathieu Desbrun, Mark Meyer, Peter Schröder, and Alan H. Barr. 2000. AnisotropicFeature-Preserving Denoising of Height Fields and Bivariate Data. In
Proc. GraphicsInterface . 145–152.Julie Digne and Carlo de Franchis. 2017. The Bilateral Filter for Point Clouds.
ImageProcessing On Line
SIAM Journal onImaging Sciences
9, 1 (2016), 116–151.Shachar Fleishman, Iddo Drori, and Daniel Cohen-Or. 2003. Bilateral Mesh Denoising.In
Proc. SIGGRAPH . 950âĂŞ953.Wendell H. Fleming and Raymond Rishel. 1960. An integral formula for total gradientvariation.
Archiv der Mathematik
11, 1 (01 Dec 1960), 218–222.Guy Gilboa. 2013. A Spectral Approach to Total Variation. In
Scale Space and VariationalMethods in Computer Vision . Springer Berlin Heidelberg, 36–47.Guy Gilboa. 2014. A Total Variation Spectral Framework for Scale and Texture Analysis.
SIAM Journal on Imaging Sciences
7, 4 (2014), 1937–1961.Guy Gilboa, Michael Möller, and Martin Burger. 2015. Nonlinear Spectral Analysisvia One-Homogeneous Functionals: Overview and Future Prospects.
Journal ofMathematical Imaging and Vision
56 (2015), 300–319.Lei He and Scott Schaefer. 2013. Mesh denoising via L minimization. ACM TOG
32, 4(July 2013), 1.Jin Huang, Tengfei Jiang, Zeyun Shi, Yiying Tong, Hujun Bao, and Mathieu Desbrun.2014. ℓ -based construction of polycube maps from complex shapes. ACM TOG
ACM TOG
22, 3 (July 2003), 943âĂŞ949.Michael Kazhdan, Jake Solomon, and Mirela Ben-Chen. 2012. Can mean-curvature flowbe modified to be non-singular?
Comput Graph Forum
31, 5 (2012), 1745–1754.Bruno Lévy. 2006. Laplace-beltrami eigenfunctions towards an algorithm that" under-stands" geometry. In
Proc. SMI . IEEE, 13–13.Peter Lindqvist. 2006. Notes on the p-Laplace equation.
Report. University ofJyvÃďskylÃď Department of Mathematics and Statistics
102 (01 2006).Hsueh-Ti Derek Liu and Alec Jacobson. 2019. Cubic Stylization.
ACM TOG
38, 6 (Nov.2019).Hsueh-Ti Derek Liu, Alec Jacobson, and Keenan Crane. 2017. A Dirac Operator forExtrinsic Shape Analysis.
Comput. Graph. Forum
36, 5 (Aug. 2017), 139–149.Simone Melzi, Jing Ren, Emanuele Rodolà, Abhishek Sharma, Peter Wonka, and MaksOvsjanikov. 2019. ZoomOut: Spectral Upsampling for Efficient Shape Correspon-dence.
ACM TOG
38, 6 (Nov. 2019).Niloy J Mitra and An Nguyen. 2003. Estimating surface normals in noisy point clouddata. In
Proc. SoCG . 322–328.Michael Moeller, Julia Diebold, Guy Gilboa, and Daniel Cremers. 2015. Learningnonlinear spectral filters for color image reconstruction. In
Proc. ICCV . 289–297.Thomas Neumann, Kiran Varanasi, Christian Theobalt, Marcus Magnor, and MarkusWacker. 2014. Compressed Manifold Modes for Mesh Processing.
Comput GraphForum
33, 5 (2014), 35–44.Stanley Osher and Leonid I. Rudin. 1990. Feature-Oriented Image Enhancement UsingShock Filters.
SIAM J. Numer. Anal.
27, 4 (Aug. 1990), 919âĂŞ–940.Maks Ovsjanikov, Mirela Ben-Chen, Justin Solomon, Adrian Butscher, and LeonidasGuibas. 2012. Functional maps: a flexible representation of maps between shapes.
ACM TOG
31, 4 (2012), 30.Fabian Prada and Michael Kazhdan. 2015. Unconditionally Stable Shock Filters forImage and Geometry Processing.
Comput Graph Forum
34, 5 (2015), 201–210.Guodong Rong, Yan Cao, and Xiaohu Guo. 2008. Spectral mesh deformation.
The VisualComputer
24, 7-9 (2008), 787–796.Leonid I. Rudin, Stanley Osher, and Emad Fatemi. 1992. Nonlinear Total VariationBased Noise Removal Algorithms.
Phys. D
60, 1-4 (Nov. 1992), 259–268.Nicholas Sharp, Yousuf Soliman, and Keenan Crane. 2019. Navigating intrinsic triangu-lations.
ACM TOG
38, 4 (2019), 55.Justin Solomon, Keenan Crane, Adrian Butscher, and Chris Wojtan. 2014. A GeneralFramework for Bilateral and Mean Shift Filtering.
CoRR abs/1405.4734 (2014).Olga Sorkine. 2005. Laplacian Mesh Processing. In
Eurographics 2005 - State of the ArtReports . 53–70.Olga Sorkine, Daniel Cohen-Or, Yaron Lipman, Marc Alexa, Christian Rössl, and Hans-Peter Seidel. 2004. Laplacian Surface Editing. In
Proc. SGP . 175âĂŞ184. Gabriele. Steidl, Joachim. Weickert, Thomas. Brox, Pavel. MrÃązek, and Martin. Welk.2004. On the Equivalence of Soft Wavelet Shrinkage, Total Variation Diffusion, TotalVariation Regularization, and SIDEs.
SIAM J. Numer. Anal.
42, 2 (2004), 686–713.Xianfang Sun, Paul L. Rosin, Ralph Martin, and Frank Langbein. 2007. Fast and EffectiveFeature-Preserving Mesh Denoising.
IEEE TVCG
13, 5 (2007), 925–938.Tolga Tasdizen, Ross Whitaker, Paul Burchard, and Stanley Osher. 2002. Geometricsurface smoothing via anisotropic diffusion of normals. In
Proc. VIS . 125–132.Gabriel Taubin. 1995. A Signal Processing Approach to Fair Surface Design. In
Proc.SIGGRAPH . ACM, 351–358.Carlo Tomasi and Roberto Manduchi. 1998. Bilateral filtering for gray and color images.In
Proc. ICCV . Narosa Publishing House, 839–846.Weihua Tong and Xue-Cheng Tai. 2016. A Variational Approach for Detecting FeatureLines on Meshes.
Journal of Computational Mathematics
34 (01 2016), 87–112.Peng-Shuai Wang, Xiao-Ming Fu, Yang Liu, Xin Tong, Shi-Lin Liu, and Baining Guo.2015. Rolling Guidance Normal Filter for Geometric Processing.
ACM TOG
34, 6(2015).Yu Wang and Justin Solomon. 2019. Chapter 2 - Intrinsic and extrinsic operators forshape analysis. In
Processing, Analyzing and Learning of Images, Shapes, and Forms:Part 2 . Handbook of Numerical Analysis, Vol. 20. Elsevier, 41 – 115.Li Xu, Cewu Lu, Yi Xu, and Jiaya Jia. 2011. Image smoothing via L gradient minimiza-tion. In Proc. SIGGRAPH Asia . ACM Press, 1.Wotao Yin and Stanley Osher. 2013. Error Forgetting of Bregman Iteration.
J. Sci.Comput.
54 (02 2013).Hao Zhang, Oliver van Kaick, and Ramsay Dyer. 2010. Spectral Mesh Processing.
Comput Graph Forum
29, 6 (2010), 1865–1894.Huayan Zhang, Chunlin Wu, Juyong Zhang, and Jiansong Deng. 2015. VariationalMesh Denoising Using Total Variation and Piecewise Constant Function Space.
IEEETransactions on Visualization and Computer Graphics
21, 7 (2015), 873–886.Qi Zhang, Xiaoyong Shen, Li Xu, and Jiaya Jia. 2014. Rolling Guidance Filter. In
Computer Vision – ECCV 2014 . Springer International Publishing.Wangyu Zhang, Bailin Deng, Juyong Zhang, Sofien Bouaziz, and Ligang Liu. 2015.Guided Mesh Normal Filtering.
Comput Graph Forum
34, 7 (2015), 23–34.Youji Zheng, Hongbo Fu, Oscar Au, and Chiew-Lan Tai. 2011. Bilateral Normal Filteringfor Mesh Denoising.
IEEE TVCG
17, 10 (2011), 1521–1530.Saishang Zhong, Zhong Xie, Weina Wang, Zheng Liu, and Ligang Liu. 2018. Meshdenoising via total variation and weighted Laplacian regularizations.
ComputerAnimation and Virtual Worlds
29, 3-4 (2018), e1827.
A MATHEMATICAL DERIVATIONSA.1 Proof of Theorem 1
For functions defined on Euclidean spaces (as opposed to surfaces,which we prove here), we refer to [Fleming and Rishel 1960].
Theorem . Let E ⊂ M be a measurable closed subset of surface M , with smooth boundary ∂ E , and let χ E be its indicator function.Then: TV ( χ E ) = length ( ∂ E ) . Proof. We split the proof by first proving TV ( χ E ) ≤ length ( ∂ E ) and then TV ( χ E ) ≥ length ( ∂ E ) . Let V be the tangent vector fieldthat attains the supremum in the weak definition of TV (see Eq. (4)),and let (cid:174) n be the unit normal vector to ∂ E . We have: TV ( χ E ) = ∫ M χ E div Vdx = | ∫ M χ E div Vdx | = | ∫ E div Vdx | We now apply the divergence theorem and Cauchy-Schwarz in-equality: | ∫ ∂ E ⟨ V , (cid:174) n ⟩ d ℓ | ≤ ∫ ∂ E |⟨ V , (cid:174) n ⟩| d ℓ ≤ ∫ ∂ E d ℓ = length ( ∂ E ) . For the other direction, we pick an arbitrary vector field U , s.t. | U | ≤ U = (cid:174) n on ∂ E . This vector field is guaranteed to existsince ∂ E is smooth. Then, applying divergence theorem: ∫ M χ E div Udx = ∫ E div Udx = ∫ ∂ E ⟨ U , (cid:174) n ⟩ d ℓ = ∫ ∂ E d ℓ = length ( ∂ E ) . ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020. onlinear Spectral Geometry Processing via the TV Transform • 199:15
Finally, from Eq. (4) we get that:length ( ∂ E ) = ∫ M χ E div Udx ≤ sup V : M→TM| V | ∞ ≤ ∫ M χ E div Vdx = TV ( χ E ) . Therefore, it follows that: TV ( χ E ) = length ( ∂ E ) . □ B ALGORITHMICSB.1 Inverse scale space
The inverse scale space method starts from the projection of theinput signal u into the kernel of the TV (yielding the mean of u )and converges in an inverse fashion to u according to the PDE: ∂ t p ( t ) = u − u ( t ) p ( t ) ∈ ∂ TV ( u ( t )) p ( ) = t → / t to account for the inverse nature of the flow, see[Burger et al. 2016, Proposition 5]). The full algorithm is summa-rized below. Other than reaching stationarity in an inverse fashion,the main difference with the forward flow is that the solution hasa piecewise-constant behavior in time, as opposed to the linear be-havior for the forward case, as seen in Section 4.3. This implies thatit is sufficient to take a single derivative of the primal variable u ( t ) to compute a spectral component ϕ t (line 7 of Algorithm 3). Algorithm 3:
Inverse scale decomposition input : signal u , max diffusion time α , no. of time steps N output : spectral representation ϕ t begin v ( ) ←− u ( ) ←− mean ( u ) for t = N do u ( t + ) Alд . ←−−−−− arg min u TV ( u ) + α ∥ u − ( u + v ( t ))∥ v ( t + ) ←− v ( t ) − ( u − u ( t + )) output ϕ t ( t ) ←− u ( t + ) − u ( t ) end end B.2 Forward vs inverse flow and the role of α We show the discretization of the forward and inverse flows andthe computation of the related time scales, to provide a better un-derstanding of the inverse relation w.r.t. time between the two andto give better insight on the role of the regularization parameter α . The inverse flow is given by Eq. (33). Its discretization yields: p ( t k + ) − p ( t k ) ∆ t k = u − u ( t k + )⇔ = ∆ t k (cid:18) u ( t k + ) − u − ∆ t k p ( t k ) (cid:19) + p ( t k + ) which – due to p ( t ) ∈ ∂ TV ( u ( t )) – is the optimality condition to:min u ∆ t k (cid:13)(cid:13)(cid:13)(cid:13) u − u − ∆ t k p ( t k ) (cid:13)(cid:13)(cid:13)(cid:13) + TV ( u ) . On the other hand, the update formula in our Algorithm 3 solves: u new PDHG ←−−−−−− arg min u α ∥ u − ( u + v ( t ))∥ + TV ( u ) . Therefore, we identify α with 1 / ∆ t k , and v ( t ) with ∆ t k p ( t k ) . Forthe latter it is important to look at the update equation for p ( t k + ) : p ( t k + ) = p ( t k ) − ∆ t k ( u − u ( t k + )) . Dividing this equation by ∆ t k + yields:1 ∆ t k + p ( t k + ) = ∆ t k + p ( t k ) − ∆ t k ∆ t k + ( u − u ( t k + )) , = ∆ t k ∆ t k + (cid:18) ∆ t k p ( t k ) (cid:19) − ∆ t k ∆ t k + ( u − u ( t k + )) . Now inserting our relation of v and p we obtain: v ( t k + ) = ∆ t k ∆ t k + v ( t k ) − ∆ t k ∆ t k + ( u − u ( t k + )) . This is the update formula in line 6 of Algorithm 3 in the case ofvariable step sizes. Finally the discretization in Eq. (33) means that: t k + = t k + ∆ t k = k (cid:213) i = ∆ t k . Thus, considering quantities on a continuous time scale we interpret: u ( t k + ) = u ( k (cid:213) i = ∆ t k ) . Interestingly, the inverse scheme is extremely forgiving in terms oftruncation errors given by the finite number of iterations in Alg. 2,and does not accumulate them as proven by Yin and Osher in [2013].Therefore, one could also consider to trade computational efficiencyagainst error by relying on the latter property of the inverse flow.
The forward flow (Eq. 8) is given by: ∂ u ∂ t ( t ) = − p ( t ) , p ( t ) ∈ ∂ TV ( u ( t )) , u ( ) = u . The implicit Euler discretization of such a flow yields: u ( t k + ) − u ( t k ) ∆ t k = − p ( t k + ) ⇒ = ∆ t k ( u ( t k + ) − u ( t k )) + p ( t k + ) which is the optimality condition to:min u ∆ t k ∥ u − u ( t k )∥ + TV ( u ) . Comparing this to our Algorithm 1, we have to identify α with ∆ t k .Again, the discretization yields: t k + = t k + ∆ t k = k (cid:213) i = ∆ t k . ACM Trans. Graph., Vol. 39, No. 6, Article 199. Publication date: December 2020.
B.3 Edge-based gradient
The gradient operator for piecewise-constant signals defined onmesh triangles can be defined, for each pair of adjacent triangles,as the jump discontinuity across the shared edge. The operator isassembled as a | E | × | F | matrix G e , with the same zero pattern asthe edge-to-triangle adjacency matrix; it contains ± D e = − T − G ⊤ e A e , where A e is a diagonal matrix of edge lengths, i.e. the area elementsfor this function space, and T is the matrix of triangle areas. B.4 Parameters and stability
We expose four parameters: the step size α , the number of spectralcomponents N (both used in Algorithm 1), and the interior stepsizes σ , τ in Algorithm 2. The only parameter which requires tuningis α , while the others can be set automatically, accordingly. Step size α . Since Algorithms 1 and 3 are unconditionally stable ,we can choose arbitrarily small or large values for α . This deter-mines the resolution of the spectral representation; to get a goodseparation of features, one should increment α according to a non-uniform sampling of the spectral domain, with higher samplingdensity within the bands having stronger spectral activity. In prac-tice, following Algorithm 3, we proceed by estimating the maximumdiffusion time by choosing the smallest value for α that makes thefidelity term negligible in the following energy:12 α ∫ M ( u ( x ) − mean ( u )) dx + TV ( u ) . For such an α , minimizing the energy above corresponds to minimiz-ing the total variation alone. This ensures that no spectral activityexists beyond the diffusion interval ( , α ) . Once the maximum possi-ble α is estimated, we rescale this value as α ←− Cα , at each iteration(with C ∈ [ , ] , typically C = . α by a differentschedule would still ensure a perfect reconstruction of the inputsignal; no details are ever lost in the analysis and synthesis process.Instead, the specific choice affects the resolution at which the spectralcomponents are extracted, and depends on the specific application,as briefly discussed in the limitations paragraph of Section 8.1. Number of spectral components N . This parameter can be deter-mined automatically, since it is equivalent to the number of rescal-ings of α required to cover the interval ( , α ) , assuming α is chosenaccording to the strategy above. Interior step sizes σ , τ . To ensure stability, these parameters mustbe chosen so as to satisfy: τσ ≤ min (∥ e ∥ ∈ E )∥ G ∥ o , where ∥ G ∥ o = sup (cid:110) ∥ Gu ∥∥ u ∥ u ∈ L (M) (cid:111) is the operator norm, andthe numerator denotes the minimum edge length on the mesh. Weestimate ∥ G ∥ o with the normest function in Matlab and set σ = τ .Âť Mesh quality.
Like most geometry processing algorithms, ourspectral decomposition is susceptible to low-quality input meshes.
Example |V| |F| N Time Rec Time
Figure 4 6079 2000 50 12.9 s -Figure 6 2590 2000 40 3.9 s -Figure 14 23308 45953 30 16.59 s 0.16 sFigure 15 44313 88622 25 31.27 s 0.46 sFigure 1 157056 314112 30 180.5 s 2.50 sFigure 20 50002 100000 20 23.69 s 0.41 sFigure 9 50002 100000 20 28.5 s 0.52 s
Table 3. Table of computational times for spectral TV. Left to right: input,number of vertices, number of faces, number of spectral components com-puted, total time for computing the spectral decomposition, total time forrecovering vertex coordinates, when needed (i.e., solution of Eq. (20) ). Tiny angles cause an increase in the gradient norm ∥ G ∥ o , lead-ing to a decrease in convergence speed of the diffusion process.To address this, we adopt theapproach of [Sharp et al. 2019],which decouples the triangula-tion used to describe the domainfrom the one used in the de-composition algorithm. In theinset figure, we show an exam-ple where degenerate triangles(visible from the bottom of thedrill, bottom row) cause distor-tion when we run the TV flowon the normal field (blue curve). With intrinsic triangulations (redcurve), we get the expected results. C RUNTIME STATISTICS
The time performance of the spectral decomposition depends on afew factors: complexity of the domain, time steps, and number ofspectral components to compute. A single iteration of Alg.(2) has aconvergence rate of k w.r.t. the primal-dual gap (see [Chambolleet al. 2010]). The number of iterations to bring the gap error belowa desired threshold however depends on the internal step sizes σ , τ (see Appendix B.4). To provide a good intuition for the actual perfor-mance we report in Table 3 the runtimes for some of the examplesshowcased in the paper. The computations were done in Matlab withan unoptimized single-threaded implementation, setting the errorthreshold of PDHG to 10 − and a maximum number of iterations to1000 and 300, for scalar and vector signals respectively.and a maximum number of iterations to1000 and 300, for scalar and vector signals respectively.