Integral Invariants for Robust Geometry_期刊版

Integral Invariants for Robust Geometry ProcessingHelmut PottmannTU WienJohannes WallnerTU GrazQi-Xing Huang Stanford UniversityYong-Liang Yang Tsinghua UniversityAbstractDifferential invariants of curves and surfaces such as curvatures and their derivatives play a central role in Geometry Processing.They are,however,sensitive to noise and minor perturbations and do not exhibit the desired multi-scale behaviour.Recently,the relation-ships between differential invariants and certain integrals over small neighborhoods have been used to define efficiently computable in-tegral invariants which have both a geometric meaning and useful stability properties.This paper considers integral invariants defined via distance functions,and the stability analysis of integral invari-ants in general.Such invariants proved useful for many tasks where the computation of shape characteristics is important.A prominent and recent example is the automatic reassembling of broken objects based on correspondences between fracture surfaces.Keywords:geometry processing,curvature,integral invariant,sta-bility,3D shape understanding1IntroductionLocal shape analysis of curves and surfaces usually employs con-cepts of elementary differential geometry like curvatures (see e.g.[do Carmo 1976;Porteous 2001]).Likewise,global shape under-standing benefits from differential geometry concepts like principal curvature lines or crest lines (see for instance [Alliez et al.2003;Hildebrandt et al.2005;Kim and Kim 2005;Ohtake et al.2004;Yokoya and Levine 1989]).However,the actual computation of curvatures for real data,given as triangle meshes or voxel grids,is a nontrivial task,because numerical differentiation is sensitive to noise.A standard method to deal with rough data is denoising and smoothing prior to numeric computation.These techniques come in two categories:Global methods which employ appropriate geo-metric flows (cf.[Bajaj and Xu 2003;Clarenz et al.2004b;Osher and Fedkiw 2002])and local ones,which use approximation by smooth surfaces (cf.[Taubin 1995;Cazals and Pouget 2003;Gold-feather and Interrante 2004;Ohtake et al.2004;Razdan and Bae 2005]).We especially want to mention the method of tensor voting (cf.[Tong and Tang 2005]).Semi-differential invariants in the sense of [Van Gool et al.1992]are a way of avoiding higher derivatives by combining reference points with first order derivatives.An alternative approach to differential geometry foremploy an exact theory of discrete analogues of ties,instead of numerically approximating the area of research,which could be called discrete etry in a narrower sense,has been investigated for a many results have been achieved –see e.g.work by and Zalgaller 1967],[Cheeger et al.1984],[Pinkall 1993],[Bobenko and Pinkall 1996],[Polthier 2002],2002],[Cohen-Steiner and Morvan 2003],[Bobenko 2005].It is possible to deal with noisy data using theories,as shown by [Hildebrandt and Polthier Klette 2006],and [Rusinkiewicz 2004].However,data appears to be neither the main strength nor the of an entirely discrete theory.Typically a discrete theory of curvatures associates them to vertices or edges or faces in a mesh,and such a discrete curvature often to appear in Computer Aided Geometric Designhas an interpretation,e.g.as ‘curvature concentrated in a vertex’,or ‘integral of curvature over a face’.This measure-theoretic in-terpretation of curvatures should not be confused with the integral invariants of the present paper.Our approach to the numerical problems inherent in the computa-tion of higher order differential invariants of noisy geometry is the following:For given 3D data,we integrate various functions over suitable small kernel domains like balls and spheres,which yields integral invariants associated with each kernel location.These in-variants turn out to have a geometric meaning and can be used as curvature estimators.This method has been initiated by [Manay et al.2004]and [Connolly 1986]and is also the topic of [Yang et al.2006;Pottmann et al.2007].The following properties of curvature estimators are important:•Robustness with respect to noise,including discretization ar-tifacts;•Multi-scale behaviour,i.e.,adaptability to the choice of reso-lution.Integration,which is an essential ingredient in the definition of in-tegral invariants,has a smoothing effect and achieves stability and robustness without the need for preprocessing.The present paper studies robustness aspects of integral invariants,as well as integral invariants related to distance functions.Thus we establish the theoretical explanation of robustness properties en-countered in numerical experiments and geometry processing algo-rithms.For the volume descriptor and similar invariants the rela-tion to shape (i.e.,curvatures)have already been established;we here present this analysis for geometry descriptors based on dis-tance functions.1A Prior WorkThe first to introduce integral invariants were [Manay et al.2004].One example is the area invariant suitable for estimating the curva-ture of a curve c at a point p ,where c is assumed to be the boundary of a planar domain D (see Fig.1,left):Consider the circular disk B r (p )of radius r and center p ,and compute the area A r (p )of the centered in a point p .The area invariant (left)is the surface area of the intersection of the disk p +rB with the domain D ,whereas the Connolly function (right)is the perimeter of the arc (p +rS )∩D divided by r .ture yields a way to estimate curvature at scale r,because features smaller than r hardly influence the result of computation.Manay et al.show the superior performance of this and other integral in-variants on noisy data,especially for the reliable retrieval of shapes from geometric databases.A similar invariant appears in earlier work(see Fig.1,right):The angle of the circular arc∂B r(p)∩D has been used by[Connolly 1986]for molecular shape analysis.If we multiply this so-called Connolly function by the kernel radius r,we obtain the length of the circular arc,which happens to be the derivative of the area invariant with respect to the kernel radius.The extensions of area invariant and Connolly function to surfaces in three-space are straightforward.One arrives at the volume de-scriptor,whose relation to mean curvature is derived by[Hulin and Troyanov2003],and which is used by[Gelfand et al.2005]for global surface matching.Its derivative with respect to r is the sur-face area of the spherical patch∂B r(p)∩D(see[Connolly1986]). The precise relation between these integral invariants on the one hand,and the curvature of planar curves and the mean curvature on surfaces on the other hand,has been derived by[Cazals et al.2003]. Principal component analysis of the domain B r(p)∩D is done by integrating coordinate functions and their products over that do-main.Therefore also principal moments of inertia and principal directions of B r(p)∩D are integral invariants.A discussion of their relation to principal curvatures and of computational issues is given by[Pottmann et al.2007],while[Yang et al.2006]shows applications and compares this method with other ways of estimat-ing principal curvatures.We also point to[Clarenz et al.2004b; Clarenz et al.2004a],who use principal component analysis of sur-face patches for feature detection.Results of a similarflavour,without the emphasis on computability and robustness,are the formulae of J.Bertrand and V.A.Puiseux (1848)which relate Gaussian curvature to perimeter and area of geodesic disks(cf.p.127of[Strubecker1969]).1B Organization and main contributions of the present paperOur paper is organized as follows:After introducing notation and some basic facts in sections2and3,we briefly discuss integral in-variants for planar curves in Section4.Those are mainly the area invariant and its derivatives with respect to kernel radius and kernel center,respectively.The analogous but slightly more involved dis-cussion of the3D counterparts,namely the volume descriptor and its derivatives,is performed in Section5.Section6studies invari-ants for curves in surfaces and shows how to obtain the geodesic curvature by integration.Invariants based on distance functions are the topic of Section7–both for surfaces and for space curves. Section8analyzes the stability of some important integral invari-ants in the presence of noise or surface perturbations.Section9 is concerned with efficient computation and implementation issues. Applications are briefly surveyed in Section10.We conclude our paper with Section11,which contains pointers to future research. The main contributions of the present paper are(i)new facts about asymptotics expansion of integral invariants and their relations to curvature,notably those computed from distance functions(ii)a thorough theoretical stability and robustness analysis of integral in-variants,with a focus on the volume descriptor;(iii)methods for computing integral invariants,especially the octree-based approach of Section9B.2Basics2A Notation and definition of integral invariantsIn this paper we assume that a curve in R2or a surface in R3is theboundary∂D of a domain D in R n,with n=2,3,respectively(locally,this is always the case,so this is no actual restriction).Wewrite1D for the indicator function of that domain:1D(x)=1ifthe point x is contained in D,and1D(x)=0otherwise.B denotes the unit ball,and p+rB denotes the ball with radius r and center p.In R2,such a ball is actually a disk,but we use thesame notation regardless of dimension.Further,the unit circle ofR2and the unit sphere of R3are denoted by S.S is the boundary of B.The symbol p+rS denotes the sphere or circle with centerp and radius r.In many cases,integral invariants,evaluated at the boundary pointp of a domain D,have the form of one of the two following convo-lution integrals:I r(p)=Zp+rBg(x)w(p−x)d x,I r(p)=Zp+rSg(x)w(p−x)d x,(1)where g is a function associated with the domain(e.g.,its indica-tor function or the distance from the boundary),and w is a weight function(e.g.,a constant).The symbol d x has various meanings, depending on the domain of integration.E.g.in dimension n=3, when integrating over the domain p+rB,it means a volume inte-gral.In dimension n=2,when integrating over p+rS,it means an arc length integral.In the remaining cases(n=3,integral over a sphere,and n=2,integral over a disk)d x means an area integral. As an example,both the area invariant and the Connolly function (cf.Fig.1)have the general form of Equation(1):we let g(x)= 1D(x)and w(x)=1.2B Integral invariants as functions of the kernel radiusFor geometry processing applications,the multi-scale behaviour of an integral invariant defined by(1)is important,which means the change of its value,if the kernel radius r varies.This leads us to consider integral invariants as univariate functions of the kernel ra-dius.A useful piece of information on these functions is the general relationddrI r(p)=I r(p),(2) which follows immediately from the product formula for integrals: I r(p)=R rI ρ(p)dρ(see[Pottmann et al.2007]).Another im-portant topic is the asymptotic behaviour when the kernel radius r tends to zero.[Pottmann et al.2007]give a general recipe for computing thefirst few terms in the Taylor expansion of I r(p).In-variants which are of this type are the area and volume functionals (see below),geometry descriptors based on the distance function (introduced in the present paper),and quantities used in principal component analysis(cf.[Yang et al.2006;Pottmann et al.2007]).2C Integral invariants as functions of the kernel centerWhen the radius used in the definition of an integral invariant I r(p) or I r(p)is kept constant,then this invariant is a function of the point ually we are interested in the values of the invariant when the point p is situated on the surface under investigation.The relations between integral invariants and geometric characteristics of the surface(mostly curvatures),which one is interested in,are no longer valid if p is not contained in the surface.Nevertheless we are interested in the behaviour of the invariant when p leaves the surface.The main reason for this is that in ac-tual computations we are often concerned with imprecise or quan-tized data,and the point for which integral invariants are eval-uated may actually lie at some distance from the hypothetical smooth surface under consideration.In order to estimate the ef-fect of these perturbations,we compute the gradient vector ∇I r =(∂∂p 1,∂∂p 2,∂∂p 3)·R p +rB g (x )d x ,and the same for I r (p ).The magnitude of perturbation is then measured by I r (p +∆p )=I r (p )+ ∆p ,∇I r (p ) +O (2),where O (2)denotes second or-der terms and , is the scalar product of vectors.Consequently,we can give a simple estimate for the change ∆I r in the integral invariant in terms of the change in the point p :Approximately, ∆I r ∆p · ∇I r .3Facts about curves and surfacesThis material is found e.g.in the monographs by [do Carmo 1976]or [Spivak 1975].References for the facts on distance functions quoted below are [Ambrosio and Mantegazza 1998]and [Pottmann and Hofer 2003].3A Planar curvesFor every point p of a sufficiently smooth curve c we can choose a Cartesian coordinate system with p as origin,such that the x 1axis is tangent to the curve (the Frenet frame).With respect to such coordinates,the curve may be written as the graph of a functionx 2=f (x 1),with f (x 1)=αx 21+βx 31+γx 41+O (x 51).With the well known formula κ=f (1+(f )2)−3/2for the curva-ture we can relate the coefficients in the Taylor expansion with the derivatives of the curvatures;the result isx 2=κ2x 21+κ 6x 31+κ −3κ324x 41+O (x 51).(3)Here the derivatives κ and κ of the curvature are with respectto x 1.For x 1=0this is also the derivative with respect to arc length.If the curve is of lesser smoothness,this Taylor expansion terminates not with O (x 51),but earlier.Because we need it later,we record in this place the coordinate of the intersection points c +r andc −r of this curve with the circle x 21+x 22=r 2(see Fig.2):c +r =(r,κ2r 2)+(−κ28,κ 6)r 3+(−κκ 12,κ 24)r 4+O (r 5).The coordinates of c −r are found by substituting −r for r in the previous formula.3B SurfacesNext,we discuss coordinate systems for surfaces.For every point p of a sufficiently smooth surface Φthere is a coordinate system with p as origin,such that the surface can be written as the graph of a functionx 3=1(κ1x 21+κ2x 22)+R (x 1,x 2),(4)where κ1,κ2are the principal curvatures of the surface in the point p ,and the remainder term R (x 1,x 2)is bounded by |R (x 1,x 2)|≤C ·(p x 21+x 22)3,i.e.,it is of third order.Mean curvature H andGaussian curvature K are defined by H =κ1+κ2,K =κ1κ2.If the surface is the boundary of the domain D ,we assume that the positive x 3axis points to the inside of D ,so that a convex domain gets nonnegative curvatures.The distance of a point x ∈R 3from the surface can be given a sign,depending on whether x is inside D or outside:dist(x ,Φ)>0⇐⇒x ∈D.Recall that the distance,signed or not,fulfills the eikonal equation ∇dist(x ,Φ) =1.ATaylor approximation of the signed distance function is given by dist(x ,Φ)=1(κ1x 21+κ2x 22)−x 3+O (3),(5)where O (3)means a third order remainder term (we use f (x )=O (k )as an abbreviation of f (x )=O ( x k )as x →0).Note that this Taylor expansion does not contain any quadratic terms which involve x 3(cf.[Ambrosio and Mantegazza 1998;Pottmann and Hofer 2003]).3C Space curvesA space curve c (u )has several orthonormal frames associated with it.For the purposes of this paragraph,we assume that u is an arc length parameter,and a dot indicates differentiation with respect tou .Then the Frenet frame {t ,h ,b }is defined by t =˙c,˙t =κh ,and b =t ×h .Here κis the curvature.Further,˙b=−τh ,where τis the torsion of the curve.By rotation of the Frenet frame about the unit tangent vector t we get the class of frames {t ,e 1,e 2}with e 1=cos φh +sin φb and e 2=−sin φh +cos φb .If wechoose the function φsuch that ˙φ=−τ,then both ˙e 1and ˙e 2are proportional to t ,and {t ,e 1,e 2}is called a rotation-minimizing frame.Now consider the ruled surface Ψ1parametrized by g 1(u,v )=c (u )+v e 1(u ).By differentiation we see that the partial deriva-tives g 1,u ,g 1,v and therefore the tangent plane of Ψ1in the point g 1(u,v )is spanned by t (u )and e 1(u ).As there is no dependence on v ,the surface is developable.An analogous result is true for the surface Ψ2which is defined via e 2.It follows that the planes orthogonal to c are orthogonal to both Ψ1and Ψ2,and thereforedist(x ,c )2=dist(x ,Ψ1)2+dist(x ,Ψ2)2.(6)As the surface Ψ1is developable,its principal frame in the pointc (u )=g 1(u,0)is {t ,e 1,e 2},with e 2as normal vector.By Meusnier’s theorem,the principal curvatures have the values κ1=κcos (h ,e 2)=κcos(φ+π/2)and κ2=0.For the surface Ψ2,the principal frame is {−t ,e 2,e 1},and the principal curva-tures have the values κ1=κcos (h ,e 1)=κcos φand κ2=0.This information will be useful when computing distance functions.3D Curves in surfacesA curve c contained in a surface Φhas an associated Darboux frame {t ,e 2,n },where t is the unit tangent vector of c ,n is the sur-face normal vector,and e 2=n ×t .If the curve is traversedwith unit speed,then ˙t=κg e 2+κn n ,˙e 2=−κg t +τg n ,and ˙n=−κn t −τg e 2.The coefficient functions κg ,κn ,τg which occur here are the geodesic curvature,normal curvature,and geodesic torsion of the curve c w.r.t.Φ,respectively.In this place we would like to mention the famous Gauss-Bonnet formula:For a closed curve c with interior D ⊂Φ,we have the identity H c κg =2π−R D K (x )d x ,with K as the Gaussian curvature of the surface Φ.4Simple invariants for planar curves4AArea invariant and Connolly functionGiven a planar curve c which occurs as the boundary of the planar domain D ,and a point p ∈c ,[Manay et al.2004]define the area invariant as the invariant I r according to Equation (1)with g =1D and w (x )=1=const .,i.e.,A r (p ):=Zp +rB1D (x )d x .(7)pp c(a)(b) p(c)(d)Figure2:Deriving the gradient of the area invariant.(a)The point p is moved towards p+d.(b)This is equivalent to moving the curvec in the opposite direction.(c)The area difference is highlighted.(d)The apparent length L d of the chord c+r−c−r when viewed in direction d contributes to the area difference.A r is the area of the intersection D∩B r(p)of the curve’s interior and the kernel disk B r(p).The perimeter of the circular arc D∩(p+rS)defines the invariantCA r(p):=Zp+rS 1D(x)d x=ddrA r(p),(8)(see Fig.1).The differential relation follows from(2).The name “Connolly function”is given to CA r(p)/r.It has been shown by [Cazals et al.2003]that there is the Taylor expansionCA r=πr−κr2+O(r3).(9) By integrating(9),we getA r=π2r2−κ3r3+O(r4).(10)This is a more precise estimate than A r≈r3arccos(rκ/2),which was given by[Manay et al.2004],in the sense that these two ex-pressions have different Taylor expansions as r→0.It is worth noting how the behaviour of both the area invariant and the Connolly function changes if the curve under consideration is not smooth.Of course,formulas(9)and(10)are useless.Sup-pose that the curve c is still smooth,but consists of two curvature continuous pieces joined together at the point p.If left and right limit curvatures have valuesκ−andκ+,then Equation(9)is obvi-ously still valid,with the arithmetic mean of the left and right hand curvature instead ofκ.If the curve is not even smooth,but piece-wise smooth with an opening angleαdifferent from180degrees, then the perimeter of the circular arc which lies inside the curve is shortened or lengthened accordingly,and we arrive atCA r=αr−κ−+κ+2r2+O(r3),(11)A r=αr2−κ−+κ+r3+O(r4),(12)where the second equation is found by integrating thefirst one.The caseα=πandκ−=κ+=κyields the smooth case.4B The gradient of the area functionalFollowing the general discussion of Section2C,we are interested in the change of the area invariant A r(p),if the point p varies.We pick a direction d and consider A r(p+d).Figures2.(a)–(d)show the change in area inflicted by the change in p:We haveA r(p+d)−A r(p)=L d· d +O( d 2),(13) where L d is the apparent length of the curve segment c∩(p+rB) when viewed from direction d(see Fig.2).Denote the two points of intersection of the circle p+rS with the curve c by c−r,c+r,and consider the chord vectorc r(p)=c+r−c−r.(14) With the symbol c⊥for a rotation about90degrees,the apparent length Ld of the curve segment visible in Fig.2is now expressed as L d= c⊥r,dd.This leads to the following theorem:Theorem1The gradient of the area invariant A r(p)with respect to p is obtained by rotating the chord vector about90degrees.The gradient and its norm are expressed by∇p A r=c⊥r=(c+r−c−r)⊥,(15)∇A r(p) =2r−κ24r3+O(r4).(16)Proof.The discussion preceding the theorem implies that the change in area is∆A r(p)= c⊥r,d +O(2),so the area gradient equals c⊥r.For the computation of ∇A r = c+r−c−r we use the Frenet frame associated with the point p and the coordinates of c+r and c−r given by Section3A:We getc r=(2r−κ2r3/4+O(r4),O(r3))T(17) and compute the norm of the gradient: c⊥r = c =[(2r−κ2r3/4+O(r4))2+O(r6)]1/2=2r[1−κ2r2/4+O(r4)]1/2= 2r(1−κ2r2/8+O(r4)),by the binomial series. The proof of Theorem1shows a phenomenon which occurs often when computing with Taylor expansions:The x1coordinate of the vector c r is known up to third order,whereas the x2coordinate has only the general term O(r3)there.It would appear that we cannot compute the norm up to third order at all.Nevertheless,the computation shows that the third order term in the x2coordinate is irrelevant.5Simple invariants for surfaces5A The volume and surface area descriptors3D counterparts of the invariants A r(p),CA r(p)are the volume descriptor V r(p)and the surface area descriptor SA r(p)of a point p on the boundary surface of a domain D:V r(p)=Zp+rB1D(x)d x,(18) SA r(p)=Zp+rS1D(x)d x=dV r(p).(19)The relation of these quantities to mean curvature is discussed in [Hulin and Troyanov2003]and[Pottmann et al.2007]:V r=2π3r3−πH4r4+O(r5),SA r=2πr2−πHr3+O(r4).(20)A r,dd(p+rS)∩Φreference planeFigure3:The apparent area A r,d enclosed by a space curve(p+ rS)∩Φwhen seen in direction of the vector d is found as area enclosed by the planar curve which arises when projecting the given space curve onto a reference plane orthogonal to d.The area is endowed with a sign,depending on the orientation of the boundary curve.With the area vector a r,we have A r,d= a r,d/ d . The normalized sphere area descriptor SA r/r2has been introduced by Connolly[Connolly1986]for molecular shape analysis.For an application in the samefield it has been studied by[Cazals et al. 2003].Equation(20)can be used to estimate the mean curvature.We de-fine the mean curvature estimators e H r(p)and H e r(p)at scale r by deleting higher order terms in the Taylor expansions in(20):e H r(p)=8−4V r(p),H e r(p)=2−SA r(p).(21)In the limit r→0,both e H r(p)and H e r(p)tend to the actual mean curvature H(p).Like in the curve case,we would like to study the behaviour of the volume descriptor for surfaces which are only piecewise smooth and piecewise curvature-continuous.Consider a point p at a sharp edge,where two smooth surface patches meet.The open-ing angle of that edge shall beα.In the smooth case(no edge),α=π.Clearly,the Taylor series of the volume descriptor starts with2αr3/3,but the higher order terms are not so obvious.A more detailed analysis shows that the curvature of the edge relative to the adjoining surfaces is irrelevant for the r4term,and that the volumedescriptor has the form V r=2α3r3−π(H−+H+)8r4+O(r5).HereH+and H−are the mean curvatures to either side of the edge.The case of a sharp corner is more complicated.5B The gradient of the volume functionalWe are now interested in the gradient∇V r(p)of the volume de-scriptor with respect to p.We consider the surfaceΦwhich is the boundary of the domain D.As in Section4B,we choose a direc-tion d and investigate the difference V r(p+d)−V r(p).The2D counterpart of this analysis is illustrated by Fig.2:the difference of volumes is given byV r(p+d)−V r(p)=A r,d d +O( d 2),(22) where A r,d is the apparent oriented area of the surface patchΦ∩B r(p)when viewed from the direction d.This area is the same as the oriented area enclosed by the projection of the closed curve c r=Φ∩∂B r(p)onto any plane which is orthogonal to d.Let c r(u)be a parameterization of this curve c r on the sphere S2r(p). We consider its area vectora r(p)=12Ic rx×d x=12ZIc r(u)×˙c r(u)du.(23)It is well known that the apparent area A r,d can be expressed as A r,d= a r(p),dd,which leads to∆V r(p)= a r(p),∆p + O(2).Theorem2The gradient of the volume descriptor V r(p)with re-spect to the point p is the area vector a r(p),whose definition in terms of the intersection curveΦ∩(p+rS)is given by Equation (23):∇V r(p)=a r(p).(24) If n(p)is a unit normal vector of the surfaceΦin the point p,which points towards the inside of the domain D,then the area vector has the following Taylor expansion:a r(p)=πhr2−3H2−K8r4in(p)+O(r5).(25)Proof.Equation(24)follows directly from the discussion preced-ing the theorem.In order to investigate the behavior of a r for r→0,we express the intersection curve c r in the principal frame associated with the point p,and we employ cylinder coordinates (ρ,φ,x3),such that x1=ρcosφ,x2=ρsinφ.The point of the curve c r which lies in the planeφ=const according to (4)is found by intersecting the curve x3=12(κ1(ρcosφ)2+κ2(ρsinφ)2)+O(ρ3)with the circleρ2+x23=r.From the coordinates of c+(r)in Section3A we read off that this point in cylinder coordinates is given byρ(φ)=r−18κ2n(φ)r3+O(r4),x3(φ)=12κn(φ)r2+O(r3),whereκn(φ)=κ1cos2φ+κ2sin2φ.When we return to Cartesian coordinates,we get a point c r(φ)of the intersection curve,and the area vector is computedby integration:a r(p)=122πRφ=0[ρ(φ)cosφ,ρ(φ)sinφ,x3(φ)]T×ddφ[ρ(φ)cosφ,ρ(φ)sinφ,x3(φ)]T.Carrying out this integration yields Equation(25).6Invariants for curves in surfacesSection4A dealt with the area invariant for a planar curve c,which arises as the boundary of a domain D.The area invariant A r(p) means that part of D whose distance from p does not exceed r. The same question can be asked if both the domain D and its boundary curve c are contained in a smooth surfaceΦ.We defineA r(p,Φ):=Area((p+rB)∩Φ).(26) It is interesting that the Taylor expansion of this area invariant given by the following theorem does not feature surface curvatures in its first two terms.Theorem3The area invariant of a curve c contained in a smooth surfaceΦhas the Taylor expansionA r(p,Φ)=πr2−κgr3+O(r4),(27) whereκg is the signed geodesic curvature of the curve c(i.e.,the curvature of the projection of c onto the tangent plane in the point p).Proof.We use a coordinate frame with p as origin and x 3axis or-thogonal to Φ.Without loss of generality Φhas the parametriza-tion x 3=g (x 1,x 2).We use the notation g ,i =∂g∂x i and ρ=px 21+x 22.As both the x 1and x 2axes are tangent to Φ,g ,1=O (ρ),g ,2=O (ρ).The surface area differential equals d x =q 1+g 2,1+g 2,2dx 1dx 2=(1+O (ρ2))dx 1dx 2.A domain D in the surface has a corresponding domain e Din the x 1,x 2plane,which arises as projection of D .If the diameter of Dis of magnitude O (ρ2),then the area of D ,which is computed as R e Dd x ,obviously equals Re D dx 1dx 2+O (ρ4).It follows that we can compute the area invariant A r (p ,Φ)to an ac-curacy of O (ρ4),if we project the curve c under consideration into the x 1,x 2plane.The result now follows directly from Equation (10). We notice that the last paragraph of the preceding proof shows that the area of the surface patch Φ∩(p +rB )equalsA r (p )=r 2π+O (r 4),(28)(which is e.g.Equ.(21)of [Pottmann et al.2007]).Remark 1The area invariant for curves in surfaces employs those points of the surface whose distance from the point p ,measured in Euclidean space,does not exceed r .We could also ask for the points whose distance,measured inside the surface Φ,does not exceed r .As it turns out,Theorem 3remains unchanged.The reason for this is that the area of a geodesic disk (i.e.,the points at distance ≤r from p )differs from the area of a Euclidean disk only by a term of magnitude O (r 4):According to the formulae of J.Bertrand and V .A.Puiseux (see p.127of [Strubecker 1969]),that area has theexpansion GA r =πr 2−π12Kr 4+O (r 5).Here K is the Gaussian curvature of the surface.We should note that the cost of computing geodesic disks usually is too high to merit their use if all we want to compute is curvatures.7Invariants from distance functionsIn applications where we have access to a distance function,itmakes sense to employ this function or functions derived from it for the definition of integral invariants.We discuss several ways to do this.7A Invariants composed from the signed distanceLet dist(x ,Φ)be the signed distance function associated with a sur-face Φ:dist(x ,Φ)is positive if x lies outside the domain bounded by Φ,and negative if x lies inside.The absolute value of the signed distance equals the distance from Φin the ordinary sense.We employ the general method of Equation (1)to define the signed distance integral D r .Letting g (x )=dist(x ,Φ)and w (x )=1leads toD r (p )=Zp +rBdist(x ,Φ)d x .(29)The relation of these descriptors to shape characteristics is de-scribed as follows:Theorem 4The signed distance integral D r (p )defined by Equa-tion (29)has the Taylor expansion D r =4πH 15r 5+O (r 6).Here H is the mean curvature of the surface Φin the point p .Proof.We consider the principal coordinate frame associated with the point p and the quadratic Taylor polynomial of the signed dis-tance function given by Equation (5).Obviously,D r =R rB `−x 3+1(κ1x 21+κ2x 22)+O (ρ3)´d x ,with ρ=(x 21+x 22+x 23)1/2.The volume integral of O (ρ3)over a volume of size O (r 3)is boundedby O (r 6).Computation of the integrals of the functions x 3,x 21,x 22finally yields the expression for D r . Note that the descriptor D r is related to mean curvature,just as the volume descriptor,which moreover is more stable (see the discus-sion below).Figure 4illustrates the similarity between these two descriptors.The sphere integral SD r corresponding to D r by dif-ferentiation has the Taylor expansion 4πH 3r 4+O (r 5).D 0.8D 0.4D noise 0.4D noise0.8V 0.8V 0.4Figure 4:Comparison of the descriptor D r derived from the signed(un-squared)distance the volume descriptor V r .Both are related to mean curvature.From top left:D r for two different kernel radii;D r for two different kernel radii with artificial noise added;V r for two different kernel radii.It is clearly visible that all descriptors try to capture the same geometric property (in this case,mean curva-ture),and that the effect of adding noise is almost eliminated when choosing a bigger neighbourhood for computation.For the kernel ball size compared to the bunny,see Fig.5.7B Invariants composed from the squared distanceWe now use the square of the signed distance function for the defi-nition of integral invariants.We study the squared distance integralD 2,r (p )=Zp +rBdist(x ,Φ)2d x ,(30)and the corresponding squared sphere distance integral SD 2,r (p )defined by integration over p +rS .The following theorem de-scribes their relation to the curvatures of the surface Φ.。