Smooth shapes with spherical topology:Beyond traditional modeling,effi cient deformation,and interaction

D.Schmitter(),P.Garc′?a-Amorena,and M.Unser

c○The Author(s)2017.This article is published with open access at Springerlink.com

Smooth shapes with spherical topology:Beyond traditional modeling,effi cient deformation,and interaction

D.Schmitter1(),P.Garc′?a-Amorena1,and M.Unser1

c○The Author(s)2017.This article is published with open access at Springerlink.com

Existing shape models with spherical topology are typically designed either in the discrete domain using interpolating polygon meshes or in the continuous domain using smooth but non-interpolating schemes such as subdivision or NURBS.Both polygon models and subdivision methods require a large number of parameters to model smooth surfaces. NURBS need fewer parameters but have a complicated rational expression and non-uniform shifts in their formulation.We present a new method to construct deformable closed surfaces,which includes exact spheres,by combining the best of two worlds: a smooth,interpolating model with a continuously varying tangent plane and well-defined curvature at every point on the surface. Our formulation is considerably simpler than NURBS and requires fewer parameters than polygon meshes. We demonstrate the generality of our method with applications including intuitive user-interactive shape modeling, continuous surface deformation,shape morphing, reconstruction of shapes from parameterized point clouds,and fast iterative shape optimization for image segmentation.Comparisons with discrete methods and non-interpolating approaches highlight the advantages of our framework.

shape modeling;spherical topology; parametric surfaces;splines;differential geometry

1 Introduction

1.1 Background

The representation ofshapes with sphericaltopology has been an ongoing research topic in computer graphics for more than three decades.The principal reason is the massive demand for closed genus-zero surfaces in industrial, architectural, and animation design as well as in biomedical imaging.Designing spherical-topology models that are simultaneously optimal with respect to several different shape characteristics still remains a challenge. Depending on whether an application involves user interaction,shape deformation,or optimization schemes,different aspects of a model are more important than others.

In user-interactive applications,a fundamental requirement is the ability to intuitively manipulate the shape.Typically,this requirement presupposes an easy way to interact directly with the surface as well as to control shapes locally. The surface deformation should be stable: a small perturbation of the surface should result in a small change of the shape.Numerical stability is crucial too:a theoretical model must remain useful in practice.On the other hand, an application might involve shape deformation as an optimization process.For example, in real-time shape recognition,approximation,and segmentation,the fast evaluation of derivativeand integral-based quantities in iterative settings is required.Further,the smoothness of the surface and the number of parameters involved can also play important roles.Usually,it is impossible to find a model that is optimal with respect to all of these requirements.In practice,a compromise ismade favoring the most important needs of a specific application.Existing models are based on polygon meshes,NURBS,or subdivision.

1.2 Overview and contribution

This article presents the full theory of a model for constructing deformable shapes with spherical topology,along with applications.This work was first introduced in a condensed form as a SIGGRAPH Asia technical brief in 2016[1]. Our framework is based on interpolating control points,similar to polygon meshes,while at the same time providing a smooth surface,formulated in the continuous domain as for NURBS.The resulting surface allows local control,is everywhere differentiable,and has a continuously varying tangent plane at every point on the surface as well as a well-defined Gaussian curvature. A major contribution is an explicit formulation of necessary conditions for the poles of the sphere to remain closed and smooth when deforming. We illustrate the use of our method with severalapplications.(1)User-interactive shape modeling:as the basis functions are interpolating, the control points lie directly on the surface of the object,which facilitates intuitive shape modeling. The basis is also finitely supported,enabling local surface control;it allows us to model a broad range of shapes by deforming a single spherical surface patch.(2)Smooth surface reconstruction from parameterized point clouds:if the underlying spherical parameterization of the samples is known, they can be easily interpolated with our model to reconstruct a smooth surface.(3)Effi cient surface deformation:by exploiting the affi ne invariance of our model,we illustrate how a fast implementation of minimum-energy deformation algorithms can be achieved.(4)Fast iterative optimization of deformation algorithms:we show how the iterative evaluation of surface and volume integrals can be effi ciently implemented for real-time optimization, and provide an example ofa segmentation algorithm for 3D medical images.

Our construction involves a class of smooth nonrational basis functions that have uniform shifts, which leads to a considerably simpler formulation than for traditional parametric methods.A controlpoint-based structure allows us to use fewer parameters than polygon or subdivision methods to achieve smoothness.Examples of the use of our method are shown in Fig.1.

Fig.1 Smooth modeling of shapes with spherical topology.Top:continuous deformation of a sphere into a gargoyle;a wood texture has been added to the surface.Bottom:shapes consisting of a single surface patch,constructed by interactive deformation of a sphere.The interpolating structure of the model allows us to intuitively design surfaces that can adopt shapes beyond classical spherical topology.Our framework is inherently smooth,facilitating natural texturing.

2 Related work

2.1 Continuous closed surfaces

The most widely used technique to construct deformable spheres in the continuous domain is NURBS [2,3]which are a subfamily of T-splines[4].Parametric NURBS surfaces are based on polynomial B-splines and are defined by a set of control points which allow local shape control[3,5–8].The main reason for using rational NURBS instead of(non-rational)polynomial B-splines is that NURBS are able to exactly reproduce conic sections[9].Conceptually,this property is equivalent to reproducing trigonometric functions,which is a necessary requirement for constructing spheres.Several ways of constructing NURBS spheres exist,e.g.,by constructing quarter or half circles and exploiting the properties of tensor-product splines or by constructing surfaces of revolution[10].NURBS typically involve the explicit characterization ofnon-uniform knot vectors with double knots. A drawback of NURBS is their rational form,which leads to complicated expressions for related integrals and derivatives[7]. Furthermore,the NURBS formulation depends on additional weight parameters,which have no intuitive interpretation. Other constructions to approximate sphere-like surfaces based on B-splines have been studied in Refs.[11,12],whereas in Ref.[13]an exact approach using exponentialsplines is proposed. Other models use(rational)B′ezier surfaces[14],which are also related to splines[15].

2.2 Discrete closed surfaces

Popular discrete methods are based on polygon meshes[16,17].With these models it is possible to represent shapes of arbitrary topology and hence, closed surfaces with spherical topology can be easily generated.A vast literature exists on mesh optimization,processing and discretizing continuousdomain operators(e.g.,see Refs.[18,19]).Polygon models are usually interpolating the control points coinciding with the mesh vertices;this property implies that the shape is modified by points which directly lie on the boundary of the object.Related to polygon models are subdivision methods[20, 21]used to construct surfaces[22–26].These methods are characterized by refinement operations iteratively applied to a set of points leading to a continuous limit surface with a certain regularity. Hence,subdivision can be seen as a hybrid method combining the discrete and the continuousdomain approach. Although in theory they are continuous,in practice,a finite number of iterations are applied,leading to a discrete mesh (thus, we categorize subdivision as a discrete method). As opposed to polygon mesh models,subdivision methods do not necessarily have interpolating control points. Different methods based on nonstationary refinement rules have been proposed to approximate spheres using subdivision[27–29].One drawback ofpolygon and subdivision methods is that they require a large number of parameters which can be a challenge when computational speed is important(e.g.,in finite element models[30]).

2.3 Spherical parameterization

The problem of finding a parameterization for an object with spherical topology is not trivial and has been tackled in Refs.[31,32].It is linked to the problem ofordering an unorganized set ofpoints or a point cloud,which is significantly harder in 3D than in 2D.In Ref.[33],a method is presented to generate a spherical parameterization of a closed surface in the continuous domain by expressing it in a basis of spherical harmonics.A related problem is surface reconstruction from a point cloud[34,35].

2.4 Interpolation

A widely used interpolating spline in computer graphics is the Catmull–Rom spline[36].However, its nature is polynomial and hence,it cannot be used to exactly parameterize a sphere.A variant ofthe Catmull–Rom spline used in signalprocessing is the Keys cubic convolution interpolator[37]which has been generalized by Refs.[38,39]to construct a trigonometric interpolation kernel that is able to reproduce conic sections.Other variants have been presented in Refs.[40–42].An interpolating subdivision scheme was originally introduced by Deslaurier and Dubuc[43].Variants of this scheme have been proposed in Ref.[44]which have also been used to construct conic sections[45].

3 Parametric shape representation

We use bold font for vectors and plain font for scalars,e.g.,c=(cx,cy,cz). To denote partialderivatives,we use the notation?σ(u,v)/?u = σu(u,v).

Note that throughout this article we will use the terms spherical topology and closed surface to describe the same kind of objects,namely connected surfaces without holes or boundaries. Using these terms to describe equivalent objects makes sense in computer graphics because in a digital environment,even continuous-domain objects can only be represented by a discretized approximation. However,in the field of mathematical topology a more rigorous definition of these terms would be required.

4 Tensor-product surfaces

4.1 Basics

We construct parametric shapes using integer shifts ofa(non-rational)generator function?.A 3D curve r(t)that is described by the coordinate functions x(t),y(t),and z(t)with t∈R is then represented by a linear combination of integer shifts of?as

where{c[k]=(cx[k],cy[k],cz[k])}k∈Zare the 3D control points. The model(1)may be extended to construct a separable parametric tensor product surfaceσ(u,v)with u,v∈R,represented as the component-wise product(denoted by the symbol×) of two curves r1×r2,i.e.

Based on this equation,an arbitrary non-separable surface with controlpoints c[k,l]can be constructed whose expression corresponds to the last line of Eq.(2).

4.2 Generator function?

The shapes that model(2)can produce depend on the generator?.For example,if?is a B-spline, the resulting shapes are polynomial.In our case,we are interested in generating trigonometric shapes in order to be able to construct exact spheres.For this purpose,we use the piecewise exponential generator proposed by Ref.[39],which reproduces sines and cosines.It is defined as?= β?ψ,whereβis a third order exponential B-spline,ψis an appropriate smoothing kernel and,?denotes convolution.We provide the explicit expression for?in Appendix A. The relevant characteristics of?for our construction, besides its sphere-reproduction property,are that it is twice differentiable,with bounded second derivatives,and satisfies the interpolation property ?(t=k)=δk,whereδkdenotes the Kronecker delta, t∈R,and k∈Z.Our generator constitutes a partition of unity,i.e.,Pwhich is a necessary and suffi cient condition for the represented shapes to be affi ne invariant.Because this generator depends on the number M of control points used to construct a curve r,we use the notation?Minstead of?.The support of?Mis equal to 4.

4.3 Defi nitions

As a simplification to indicate the M1-periodized basis function,we write:

andφ2:= ?2M2.To denote the integer shifts of the basis functions on the normalized parameter domain we useφ1,k(t):=φ1(M1t?k)andφ2,k(t):= φ2(M2t?k).

5 Spherical parameterization

5.1 The deformable sphere

In this section,we outline our proposed construction ofthe deformable sphere.Without loss ofgenerality, we parameterize its surface as

with u,v∈[0,1].In Ref.[39],it has been shown that:

The periodization ofφ1as defined in Eq.(3)allows us to express the u-dependent 1-periodic trigonometric functions in Eq.(4)using a finite sum and M1∈Z control points.The v-dependent trigonometric functions in Eq.(4)are not periodic and are expressed as

Because the support of?Mis equal to 4,for v∈[0,1],we have?M2(v?l)=0 if l?[?1,...,M2+ 1],which explains the limits of the sum in Eq.(6). Following the construction given in Eq.(2),we finally parameterize the sphere as

where the control points of the surface are given by its samples:

Note that M1and M2are the numbers of control points in the u-and v-directions.Hence,this representation allows us to construct a perfect sphere with any numbers M1,M2of control points.The only condition for the integer shifts of?Mto form a stable basis,i.e.,to guarantee a stable implementation,is M ≥ 3[39].A reconstructed sphere with interpolatory control points is shown in Fig.2.

5.2 Smoothness conditions at the poles

Fig.2 Reconstructed sphere with interpolatory controlpoints shown in green.The parametric directions are indicated by the blue and red arrows.

Since? ∈ C1,continuity is guaranteed nearly everywhere on the surface as long as the control points do not overlap.However,for the deformed sphere,smoothness is not guaranteed at the poles unless we take appropriate measures.In Ref.[11],it is shown that continuity at the poles is ensured ifthe deformable sphere is constructed with continuously varying tangent planes at these points.This condition is expressed mathematically as

for the South pole,where T1,N,T2,N,T1,S,and T2,Sare vector parameters that can be freely chosen.In Appendix B,we show that both sides of Eq.(9)can be simplified independently and we end up with the condition:

Similarly,Eq.(10)simplifies to

The tangent plane at the poles is then spanned by the vectors T1,N,T2,Nand T1,S,T2,S.Figure 3 illustrates the effect of imposing the smoothness conditions at the poles.

5.3 Interpolation conditions at the poles

The sphere needs to remain closed when deforming in order to maintain spherical topology.Again,special attention needs to be paid to the poles:allthe circles of longitude of the original sphere should originate and end at the poles of the surface.In accordancewith the parameterization in Eq.(4),this condition is expressed as

Fig.3 Closed and smooth deformable sphere.Left:if no smoothness conditions are imposed,the surface is non-diff erentiable at the poles. Center:if no pole-interpolation conditions are imposed,the surface looses its spherical topology when deforming.Right:a closed and deformed sphere is shown with smoothly varying tangent planes at the poles.

In Appendix C,we show that condition(13) translates directly into

?k∈ [0,...,M1?1]. In Fig.3,we compare a deformed sphere with and without imposing the closeness conditions at the poles.

5.4 Main result

We combine allof the above considerations together in order to state the main result of this article.A locally and smoothly deformable sphere is expressed by the parameterization in Eq.(7)subject to the smoothness conditions(11)and(12)and the closeness condition(14).

5.5 Useful properties ofσin practice

Our deformable sphereσis affi ne invariant.Hence, its construction is independent of location and orientation,i.e.

where A is a 3×3 matrix and b a constant vector in 3D.

Furthermore,sinceφis twice differentiable,the surface has everywhere a well-defined tangent plane and Gaussian curvature.This property,for instance, allows us to compute the normalvector at any point on the surface,an important requirement to render a textured surface.

6 Results and applications

6.1 Interactive modeling

It is crucial in interactive shape modeling that the modeling process is intuitive.Standard modeling applications allow a user to modify a shape by dragging its control points with the mouse in order to displace them.If the control points lie directly on the surface of the shape,the modeling task is significantly simplified. This is the case for polygon models,but then the underlying shape is not smooth.On the other hand,NURBS allow for the construction of smooth shapes,but the control points do not interpolate the shape.This makes the modeling task less intuitive. Local shape control is diffi cult as the surface becomes more complex because it is no longer clear which part ofthe surface is affected by a specific control point.Our proposed construction solves this problem since?Msatisfies the interpolation condition and is also smooth. Hence,even if the modeled surface is of great complexity,the modeling process remains intuitive and simple since the control points always lie on the boundary of the shape.Furthermore,thanks to the compact support of?M,local shape control is guaranteed.Figure 4 illustrates the interactive shape modeling process.

Fig.4 Interactive modeling.Left:the region(yellow)aff ected by moving a single control point(blue)is shown;it corresponds to a patch of size 4×4 due to the support of the generator?M.Right:a brain(green)modeled using our interpolatory construction(bottom) and compared to the process where a non-interpolatory basis function is used(top)similar to NURBS.The coordinate system indicates a control point about to be interactively displaced in 3D space.Top right:it is unclear which region ofthe surface is controlled by a certain control point.The two poles are indicated in the fi gure to show the importance ofthe smoothness property in practice.Center:a smooth brain model rendered based on the modeling process illustrated on the right.

6.1.1 Intuitive user interaction

Our proposed framework can be exploited to make user-interactive shape modeling more intuitive and compelling.It is ideally suited for implementation in interactive shape modeling software,where the user modifies the shape by displacing the control points with the mouse.With relatively few control points,complex structures are easily constructed and modified. Figure 5 shows examples of the use of our framework in an interactive modeling environment. Final renderings,where texture is added to a shape,are achieved without discretization artifacts,independently of the number of control points chosen,since the underlying structure is smooth(see Fig.1,bottom row).In Fig.6,the effect of constraining the poles to be smooth is illustrated when performing interactive modeling.

6.2 Shape interpolation

When dealing with a parameterized point cloud whose points correspond to the samples of a surface with spherical topology,our formulation allows for an immediate reconstruction of the smooth shape. Several algorithms have been proposed to obtain such a parameterization(see Section 2). In this case,for each point p∈R3ofthe point cloud,a pair (uk,vl)of coordinates is assigned in the parameter domain and we can establish the relationσ(uk,vl)= pk,l=c[k,l].For fixed numbers ofpoints,M1,M2,in the u-and v-directions,the parametric coordinates for the normalized parametric domain,i.e.,u,v∈[0,1],are given by uk=k/M1and vl=l/M2.The resulting continuously defined surfaceσ(u,v)is immediately reconstructed since it is fully specified by its control points subject to the smoothness and pole-interpolation conditions described above.An example is shown in Fig.7.

6.2.1 Smooth modeling at arbitrary resolution

Fig.5 Implementation of the framework in a shape modeling environment.Diff erent shapes are interactively designed starting from a sphere (from left to right).The interpolatory control points allow us to easily model surfaces that can adopt shapes beyond traditional spherical topology,such as the mug,rocket,or bullet.The last two rows show shapes where only the closeness condition has been imposed in order to allow for the construction of sharp kinks.

Fig.6 Poles with continuously varying tangent plane.The eff ect of imposing the smoothness condition on the poles in interactive shape modeling is illustrated.Left:smooth pole.Right:sharp discontinuity at the pole resulting in a singularity.

Fig.7 Interpolation of a parameterized point cloud. The dinosaurs(middle and right)are smooth reconstructions obtained by interpolating the point cloud on the left.Our surface construction is affi ne invariant and hence,rotating the shape is simply performed by rotating the point cloud.

Because our construction ofσis inherently smooth, even with few control points,the tangent plane and Gaussian curvature are everywhere well-defined. This advantageous property allows construction of textured models with few parameters,which is usefulfor example in applications involving real-time rendering.As an example,we have parameterized the point cloud of the Gargoyle model using the algorithm described by Ref.[32],which allows us to reconstruct a smooth surface by interpolating the points.Additionally,we have subsampled the point cloud at different resolutions to obtain an approximation ofthe Gargoyle with varying levels of accuracy.Figure 8 illustrates the result and makes a comparison with a model based on polygons.

6.2.2 Compression

Related to the previous example is the problem of shape compression.Typically,the fewer coeffi cients are used to compress a smooth shape,the more discontinuous its representation becomes, which influences its texturing and rendering.The advantage of our modelis that smoothness is always preserved,even with few coeffi cients,as shown in Fig.8.

6.3 Effi cient shape deformation

Fig.8 Interpolation of shapes with spherical topology:smooth Gargoyle reconstructions at diff erent resolutions.The same number of control points is used in both directions of the parameter domain, i.e.,M=M1=M2.Top:results obtained with our construction. Bottom:a(linear)polygon reconstruction method is applied.Note that with our approach the smoothness of the model does not depend on the number of parameters.

An advantage of using continuous-domain models based on control points is that the shapes are described by a finite number of control points, whereas the corresponding coordinate functions x,y, and z live in an infinite dimensional space;this allows us to describe a shape deformation process in the continuous domain just by displacing the control points. In the following,we provide two examples that illustrate how the minimization of distance criteria in the continuous domain can be effi ciently formulated as conditions on the control points. Other deformation criteria which can be minimized in a similar way have been studied in Refs.[46,47].

6.3.1 Minimum-energy deformation

We illustrate two deformation processes which correspond to minimum-energy deformation in L2([0,1]2,R3). Both processes are formulated entirely with respect to the control points.Thereby, we can parameterize the path,which describes the deformation in the space that contains all parametric shapes.An immediate application is the construction of interpolated or extrapolated shapes, where the terms interpolated and extrapolated refer to a shape lying on the path in some shape space. Typically,such a shape space is described by a metric that provides a notion of distance between two points that lie in the space.Hence,in a given shape space,a shape is treated as a single point. Here we are interested in describing the deformation such that a minimum amount of energy is required in order to deform one shape into another.This translates into describing the deformation as the shortest path between two points in the shape space according to its underlying metric.

TheHilbertplaneasashapespace.Given two surfacesσ1andσ2living in the Hilbert plane L2([0,1]2,R3),the shortest path connecting them can be parameterized by the intermediate surfaceσ that minimizes:

for a givenτ∈[0,1].We see immediately that,for τ=0,the minimizer isσ=σ2,whereas forτ=1 it isσ=σ1.For values ofτ∈R[0,1],the path F describes extrapolated shapes,i.e.,shapes that do not lie between the two surfacesσ1andσ2.The L2-norm in Eq.(15)is induced by the L2-inner product:

Using the property that our parameterization is affi ne invariant,it is easy to show that the solution of min F(τ,σ)is given by

where C,C1,C2are the matrices which contain all the control points of the corresponding surfaces. As an example that illustrates the deformation process and also the effect of imposing the closenesscondition on the poles,we have deformed a disk into a sphere. Figure 9 illustrates this process and compares it to the case where no poleinterpolation conditions are imposed.Figure 1 shows the deformation of a sphere into a Gargoyle.

TheHilbertsphereasashapespace.Every parametric shape can be projected onto the unit Hilbert sphere by normalizing it such that‖σ‖L2=

Fig.9 Minimum-energy deformation in the Hilbert plane.Top: a disk is deformed into a sphere through Eq.(17).Bottom:the same process,but without imposing the pole-interpolation conditions in Eq.(14).In this case,the surface does not remain closed when deforming and Eq.(17)describes the deformation between a circle and a sphere.

whereθ=cos?1(〈σ1,σ2〉),Γ(0)=σ1,andΓ(1)= σ2. Again, if τ∈[0,1], Eq. (18) describes interpolated shapes,whereas forτ∈R[0,1],τ describes extrapolated shapes.As in the previous example, we exploit the affi ne invariance of our parameterization in order to describe the deformation as a function of the controlpoints.The interpolating control points are given by

An example invoking this deformation is shown in Fig.10.

Morphing. Computing morphs between two or several shapes is similar to computing the deformation between shapes.The difference is that the deformation is expressed as a parameterized weighted linear combination of two shapes,whereas a morph corresponds to a particular instance of the parameterized function.Concretely,if Eq.(17)or Eq.(19)is evaluated for a specific value ofτ,we obtain a morph betweenσ1andσ2. Examples of such smooth morph constructions are shown in Fig.10,which correspond to morphed point clouds similar to the ones shown in Fig.11.

Fig.10 Minimum-energy deformation on the Hilbert sphere.Top: sphere and Venus.Bottom:Stanford Bunny and Gargoyle.Bothfrom left to right.

Fig.11 Infl uence of the parameterization on the deformation. The point clouds(with M=M1=M2=270)that define the control points of the dinosaur and the Gargoyle are parameterized and the locations of the poles are indicated with red arrows.The dinosaur(top left)is deformed into the Gargoyle(top right).The two intermediate shapes in the top row illustrate the deformation. Bottom:the poles on the sphere can be placed at different locations. For instance,if a diff erent parameterization of the dinosaur is chosen such that the the North pole cNand South pole cSare exchanged, then the deformation process is diff erent.

Parameterization. An important aspect to consider when using our model is that the parameterization which describes the shape is not unique.This is natural in the case of surfaces with spherical topology and originates from the fact that there exist an infinite number of ways to place the two poles(with the constraint that they must be opposite to each other).However, considering Eqs.(17)and(19),we see that there is a unique correspondence between two given spherical parameterizations,which implies that given two surfaces,σ1andσ2,each control point c1[k0,l0]is transformed into c2[k0,l0].If a different parameterization is chosen for at least one ofthe two surfaces(i.e.,ifthe poles are placed differently),then the resulting deformation willinevitably be different. This is illustrated in Fig.11.Insights into finding an optimalcorrespondence between shapes can be found in Refs.[48]and[49].

6.4 Fast computation of surface and volume integrals

In certain applications that require iterative optimization,it is necessary to rapidly compute surface or volume integrals effi ciently.An example is the deformation of a surface guided by optimizing an energy functional in real time.

6.4.1 Flux across surface

We illustrate how a flux E across a surface S, parameterized byσ(u,v),may be computed rapidly and effi ciently.Given a vector field f,one way of expressing the flux E is by

where d S represents the vector diff erential of the surface area,∧ denotes the wedge product, anddiv f(τ,y,z)d τis the preintegrated divergence of the vector field f in the x-direction.Typically,f does not depend on the surface and hence,gxcan be precomputed and stored in a look-up table to significantly speed up the computation.We derive Eq.(20)in Appendix D. The use of pre-integrated functions is only possible because we define the surfaceσin the continuous domain.Next,the flux E can be effi ciently optimized by computing the gradient of E with respect to the control points using a gradient-based iterative method.An explicit expression of the gradient can be obtained easily,and hence implemented in an exact way.

Example.We illustrate the above computation by segmenting the surface of a human brain in a 3D MRI image.We first compute an edge map of the 3D image using a standard surface extraction algorithm[50]and construct an energy functional E that depends on the gradient of the edgemap. Hence,in Eq.(20),the gradient becomes f.By minimizing Eq.(20),σdeforms iteratively in order to approximate the edge map,as shown in Fig.12. The result can easily be manually adjusted by a clinician(see Fig.4),which is an additional advantage of our algorithm compared to existing methods.

6.4.2 Exact volume computation

For ki∈[0,...,M1?1]and li∈[?1,...,M2+1]the volume enclosed by the surfaceσis computed by

where cx,cy,and czare the x,y,and z coordinates of the control points ofσand

Fig.12 Brain segmentation in a 3D medical MRI image.The red surface is a rendered edge map extracted from medical data. An ellipsoidal surface is initialized inside the brain surface(left)and evolves by iteratively minimizing(20)(from left to right).The fi nal result(right)corresponds to a smooth and continuous closed surface shape.

Sinceαdoes not depend on the control points c,it can be precomputed and stored in a look-up table in order to quickly evaluate the volume in interactive optimization schemes.Furthermore,becauseφand its derivativeφ′have compact support,the number of non-zero elements in Sum(21)is small,which additionally simplifies the computation.We derive the formula for the volume in Appendix E.The integrals in Eq.(22)can be further simplified and exactly evaluated using techniques from spline theory similar to the approaches in Refs.[51,52].

7 Implementation

In this section,we describe some important details regarding the implementation.

7.1 Choice of free parameters

7.1.1 Exact sphere

The orientation of the sphere is given by Eq.(4) and therefore,the coordinates of the North pole are cN=(0,0,1),and for the South pole,cS=(0,0,?1). Since by construction,the vectors T1,Nand T2,Nspan the tangent plane at the North pole of the sphere,a naturalchoice is to set T1,N=(1,0,1)?cN=(1,0,0) and T2,N=(0,1,1)?cN=(0,1,0). With the same approach we also obtain T1,S=(1,0,0)and T2,S=(0,1,0)for the South tangent plane.Note that because our construction is affi ne invariant,for a sphere with a diff erent size or orientation,the new coordinates are found by applying the corresponding affi ne transformation to the existing control points.

7.1.2 Arbitrary shape with spherical topology

A simple method to estimate the tangent plane is to compute the plane that best approximates the points lying on the first circle of latitude next to the North or South pole.Any two vectors spanning this plane can be chosen as T1,N,T2,Nand T1,S,T2,S.

7.2 Discretization of basis functions

Because our construction is formulated in the continuous domain,the shape representation can be discretized with arbitrary precision in order to implement it.An effi cient way is to discretize the interpolator?rather than the surface,which becomes highly beneficial,for example,in interactive applications where the shapes to be constructed are not known beforehand.By discretizing?prior to surface construction,the samples of the interpolator can be stored in a look-up table to speed up the surface reconstruction.Thus,the sampling rate is freely chosen.In Figs.13 and 14,we show the effects of different sampling rates.A sampling rate equal to one corresponds to a polygon model(i.e.,linear interpolation between points),which means only the blue sample in Fig.13 is non-zero.The higher the sampling rate,the closer the approximation of the continuous domain model.Its effect on surface reconstruction is shown in Fig.14.In practice,if a large number of control points is used,one can already obtain satisfactory smoothness ofthe surface with a low sampling rate,whereas with a small number of control points,the sampling rate must be higher to obtain a smooth surface.

8 Discussion and future work

8.1 Comparison with NURBS

Fig.13 Sampling of the interpolator?.Because?is formulated in the continuous domain,it can be discretized with arbitrary precision. If only one sample is considered(blue sample),then the result corresponds to linear interpolation,which is equivalent to a polygon model.The samples denoted by orange circles correspond to a lower sampling rate than the green samples.

Fig.14 Eff ects of diff erent sampling rates,increased from left to right.The surfaces are constructed with M=M1=M2.The red wireframe corresponds to a lower number of control points(M=20) used to reconstruct the bunny than the blue wireframe(M=40).

Our formulation has several advantages compared to an approach using NURBS.With the NURBS formulation,a sphere can only be represented using multiple surface patches.The NURBS formulation requires not only more control points to represent a sphere,but also more parameters in total,due to the weights used in that formulation.Further, because the basis functions of NURBS are rational, the computation of derivatives and integrals results in complicated expressions.This can become a problem when integral-dependent quantities need to be computed,such as in the evaluation of surface or volume integrals in optimization schemes. Also,the optimization itself must be carried out simultaneously with respect to the control points and the weight parameters,which introduces additional complexity.For interactive shape design, the interpolation property of our framework makes the modeling task more intuitive.Complex shapes that require more detail,and hence,more control points,are especially modeled more easily with our solution(see Fig.4).Furthermore,interpolating parameterized point clouds with spherical topology is diffi cult with NURBS due to their noninterpolatory nature;it involves complex NURBS approximation techniques or inverse filtering,which is not straightforward because of the smoothness conditions at the poles.The only NURBS that are interpolatory are zero and first degree NURBS, which are non-smooth.

Moreover,our formulation allows for a shape representation using only integer shifts.NURBS usually have non-uniform shifts.The advantage of considering integer shifts is that it allows convolution and filtering techniques as well as frequency domain calculus.This can be usefulwhen performing surface resampling,projections onto other spline spaces,and evaluation ofinner products,for example to compute L2-distances between surfaces.It also allows for a simpler formulation of the surface by specifying control points instead of non-uniform knot vectors including double knots.

8.2 Comparison with polygon models

Polygon models are inherently interpolatory schemes because the controlpoints coincide with the vertices of the mesh.Similar to subdivision schemes,these models require more parameters than our model in order to achieve a higher degree of smoothness (see Fig.8).Geometric operators and quantities, such as tangent planes,normals,curvatures,or the Laplacian have to be approximated by polygon mesh processing techniques. The same holds true for integral and derivative-based quantities.However, polygon meshes do not require an underlying parameterization of the model.

8.3 Comparison with the Catmull–Rom interpolator

The Catmull–Rom[36]or Keys interpolator[37] are interpolating and smooth.Because they are polynomial it is not possible to construct an exact sphere with these functions.However,construction of a model with spherical topology(which excludes exact spheres and ellipsoids)is possible by replacing ? in our framework with the Catmull–Rom spline.Because its support is the same as for?, our formulation for the smoothness and interpolation conditions at the poles can easily be adapted to the purely polynomial case.

8.4 Extending the framework to other topologies

Our concept can be extended to surfaces with other topologies(e.g.,cylindrical or rectangular)in order to create a unifying framework for smooth shape modeling with interpolatory control points. These topologies do not require special attention to poles and are easier to parameterize using tensor products and a suitable interpolator.One way to parameterize the rectangle is with the polynomial Keys interpolator[37],whereas the cylinder may be parameterized using? for the trigonometric part(i.e.,the circles in one direction)and the Keys interpolator for the linear part(i.e.,along the axis).Another example is the torus which is easily parameterized using?since it is periodic in u and v. In Fig.15,examples of these topologies are shown as well as how they can be smoothly deformed byexploiting the interpolation property in interactive settings.

Fig.15 Smooth modeling of diff erent topologies with interpolatory control points. Top:idealized shapes that define the topologies. The red points indicate the interpolatory control points.Bottom: a smooth deformation of these shapes.

9 Conclusions

The standard method for smooth,parametric shape modeling in industry is NURBS.In this paper,we presented an alternative method to model smooth shapes with spherical topology.The fundamental difference with the existing standard is that our basis functions are interpolatory and non-rational and we only use uniform shifts.Our formulation is simpler than NURBS and thus has several advantages in practical applications, including immediate reconstruction of smooth surfaces by interpolating parameterized point clouds,more intuitive shape modeling,and simplified formulation ofoptimization algorithms that involve integral-and derivativedependent quantities.Our framework is extensible to a richer family of topologies.A video illustrating the use of our framework in practice is available at http://bigwww.epfl.ch/demo/siggraph2016/.

Appendix AExplicit expression for?

The generator is given by?=β?ψ,where?denotes continuous convolution.The functionβis a thirdorder exponential B-spline defined by

are constants that depend only on M.The explicit expression for?=?Mis given by Eq.(23).Note that?is non-rationalwith respect to its parameter.

Appendix BSmoothness conditions at poles

The left-hand-side of Eq.(9)is developed as

where we have used the fact that?′is odd(?is even). The right-hand-side of Eq.(9)may be expressed as

By equating Eqs.(24)and(25)and identifying the coeffi cients,we obtain Eqs.(11)and(12).

Appendix CInterpolation conditions atpoles

At the North pole,we have:

Since?Msatisfies the interpolation condition,the term that depends on l is always zero unless l=0??2M2(l=0)=1.Hence,Eq.(26)simplifies to

Because the integer shifts of?M1build a basis[39] and?M1satisfies the partition of unity property, Eq.(27)only holds if c[k,0]=C,with C a constant vector for all k.Thus,

A similar derivation leads to the interpolation condition at the South pole.

Appendix DFlux across surface

We denote by n the normalvector to the surface and make use of the divergence theorem to compute:

where gx,gy,gzare the pre-integrated functions along directions x,y,or z.The wedge operator is defined as

and is explicitly computed using?σ/?u=(?x/?u,?y/?u,?z/?u)and?σ/?v=(?x/?v,?y/?v,?z/?v) with

Appendix EVolume computation

By the divergence theorem,the volume of a parametric surface is given by

By applying the same simplifications to compute the wedge operator(28)as in Appendix D,and using

Because only the basis functions depend on u and v,the integral-dependent terms can be isolated and precomputed to obtain Eq.(21).

Acknowledgements

This work was funded by the Swiss National Science Foundation under Grant 200020-162343.We are grateful to Zsuzsanna P¨usp¨oki for help with the figures and to Irina Radu for help with the video. We also appreciate the interesting discussions on the subject that we had with Masih Nilchian and Emrah Bostan.We thank Mike McCann for proof-reading the manuscript.

[1]Schmitter,D.; Garc′?a-Amorena,P.; Unser,M. Smoothly deformable spheres:Modeling,deformation, and interaction.In:Proceedings of the SIGGRAPH ASIA 2016 Technical Briefs,Article No.2,2016.

[2]Piegl,L.On NURBS:A survey.IEEE Computer Graphics and Applications Vol.11,No.1,55–71,1991.

[3]Piegl,L.;Tiller,W.The NURBS Book,2nd edn. Springer Berlin Heidelberg,2010.

[4]Sederberg,T.W.;Zheng,J.;Bakenov,A.;Nasri, A.T splines and T NURCCS.ACM Transactions on Graphics Vol.22,No.3,477–484,2003.

[5]Lyche,T.Discrete cubic spline interpolation.BIT Numerical Mathematics Vol.16,No.3,281–290,1976.

[6]Lyche,T.Local spline approximation methods. Journal of Approximation Theory Vol.15,No.4,294–325,1975.

[7]Manni,C.;Pelosi,F.;Sampoli,L.Generalized B splines as a tool in isogeometric analysis.Computer Methods in Applied Mechanics and Engineering Vol. 200,Nos.5–8,867–881,2011.

[8]Schumaker,L.L.Spline Functions:Basic Theory.J. Wiley&Sons,1981.

[9]Farin,G.E.NURBS:From Projective Geometry to Practical Use,2nd edn.Natick,MA,USA:A.K. Peters,Ltd.,1999.

[10]Rogers,D.F.Preface.In: An Introduction to fNURBSg,The Morgan Kaufmann Series in Computer Graphics.Rogers,D.F.Ed.San Francisco:Morgan Kaufmann,xv–xvii,2001.

[11]Dierckx,P.Algorithms for smoothing data on the sphere with tensor product splines.Computing Vol. 32,No.4,319–342,1984.

[12]Gmelig Meyling,R.H.J.;Pfluger,P.R.B spline approximation of a closed surface.IMA Journal of Numerical Analysis Vol.7,No.1,73–96,1987.

[13]Delgado Gonzalo,R.;Chenouard,N.;Unser,M.Spline based deforming ellipsoids for interactive 3D bioimage segmentation.IEEE Transactions on Image Processing Vol.22,No.10,3926–3940,2013.

[14]Prautzsch,H.;Boehm,W.;Paluszny,M.B′ezier and B Spline Techniques.Springer Verlag New York,2002.

[15]Romani,L.;Sabin,M.A.The conversion matrix between uniform B spline and B′ezier representations. Computer Aided Geometric Design Vol.21,No.6,549–560,2004.

[16]Botsch,M.;Kobbelt,L.;Pauly,M.;Alliez,P.;L′evy, B.Polygon Mesh Processing.CRC Press,2010.

[17]Botsch,M.;Steinberg,S.;Bischoff,S.;Kobbelt,L. OpenMesh–A generic and effi cient polygon mesh data structure.In:Proceedings ofthe OpenSG Symposium, 2002.

[18]Sorkine,O.Diff erential representations for mesh processing.Computer Graphics Forum Vol.25,No. 4,789–807,2006.

[19]Sorkine,O.;Cohen-Or,D.;Lipman,Y.;Alexa,M.; R¨ossl,C.;Seidel,H.P.Laplacian surface editing.In: Proceedings of the Eurographics/ACM SIGGRAPH Symposium on Geometry Processing,175–184,2004. [20]Dyn,N.;Farkhi,E.Spline subdivision schemes for compact sets.A survey.Serdica Mathematical Journal Vol.28,No.4,349–360,2002.

[21]Lounsbery, M.; DeRose, T.D.; Warren, J. Multiresolution analysis for surfaces of arbitrary topological type.ACM Transactions on Graphics Vol. 16,No.1,34–73,1997.

[22]Catmull,E.;Clark,J.Recursively generated B spline surfaces on arbitrary topological meshes.In:Seminal Graphics:Poineering Eff orts that Shaped the Field. New York:ACM,183–188,1998.

[23]DeRose,T.;Kass,M.;Truong,T.Subdivision surfaces in character animation.In:Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques,85–94,1998.

[24]Doo,D.;Sabin,M.Behaviour of recursive division surfaces near extraordinary points.Computer-Aided Design Vol.10,No.6,356–360,1978.

[25]Loop,C.;Schaefer,S.Approximating Catmull–Clark subdivision surfaces with bicubic patches.ACM Transactions on Graphics Vol.27,No.1,Article No. 8,2008.

[26]Stam,J.;Loop,C.T.Quad/triangle subdivision. Computer Graphics Forum Vol.22,No.1,79–85, 2003.

[27]Charina,M.;Conti,C.;Romani,L.Reproduction of exponential polynomials by multivariate non stationary subdivision schemes with a general dilation matrix.Numerische Mathematik Vol.127,No.2,223–254,2014.

[28]Novara,P.;Romani,L.Building blocks for designing arbitrarily smooth subdivision schemes with conic precision.Journal of Computational and Applied Mathematics Vol.279,67–79,2015.

[29]Warren,J.;Weime,H.Subdivision Methods for Geometric Design:A Constructive Approach.San Francisco:Morgan Kaufmann,2001.

[30]Hughes,T.J.R.The Finite Element Method:Linear Static and Dynamic Finite Element Analysis.Courier Corporation,2012.

[31]Asirvatham,A.;Praun,E.;Hoppe,H.Consistent spherical parameterization. In: Computational Science–ICCS 2005.Sunderam,V.S.;van Albada,G. D.;Sloot,P.M.;Dongarra,J.Eds.Springer Berlin Heidelberg,265–272,2005.

[32]Praun,E.;Hoppe,H.Spherical parametrization and remeshing.ACM Transactions on Graphics Vol.22, No.3,340–349,2003.

[33]Brechb¨uhler, Ch.; Gerig, G.; K¨ubler, O. Parametrization of closed surfaces for 3D shape description.Computer Vision and Image Understanding Vol.61,No.2,154–170,1995.

[34]Berger,M.;Tagliasacchi,A.;Seversky,L.;Alliez, P.;Levine,J.;Sharf,A.;Silva,C.State of the art in surface reconstruction from point clouds.In: Eurographics 2014–State of the Art Reports.Lefebvre, S.;Spagnuolo,M.Eds.The Eurographics Association, 161–185,2014.

[35]Kazhdan,M.;Bolitho,M.;Hoppe,H.Poisson surface reconstruction.In:Proceedings of the 4th Eurographics Symposium on Geometry Processing, 61–70,2006.

[36]Catmull,E.;Rom,R.A class of local interpolating splines.Computer Aided Geometric Design Vol.74, 317–326,1974.

[37]Keys,R.Cubic convolution interpolation for digital image processing.IEEE Transactions on Acoustics, Speech,and Signal Processing Vol.29,No.6,1153–1160,1981.

[38]Schmitter,D.;Delgado-Gonzalo,R.;Unser,M.A family of smooth and interpolatory basis functions for parametric curve and surface representation.Applied Mathematics and Computation Vol.272,No.1,53–63, 2016.

[39]Schmitter,D.;Delgado-Gonzalo,R.;Unser,M. Trigonometric interpolation kernel to construct deformable shapes for user interactive applications. IEEE Signal Processing Letters Vol.22,No.11,2097–2101,2015.

[40]Beccari,C.V.;Casciola,G.;Roman,L.Construction and characterization ofnon uniform localinterpolating polynomial splines.Journal of Computational Applied Mathematics Vol.240,5–19,2013.

[41]Schmitter,D.;Gaudet-Blavignac,C.;Piccini,D.; Unser,M.New parametric 3D snake for medical segmentation of structures with cylindrical topology. In:Proceedings of the IEEE International Conference on Image Processing,276–280,2015.

[42]Zhang,R.J.Uniform interpolation curves and surfaces based on a family of symmetric splines.Computer Aided Geometric Design Vol.30,No.9,844–860,2013.

[43]Deslauriers,G.; Dubuc,S.Symmetric iterative interpolation processes. In: Constructive Approximation.DeVore,R.A.;Saff,E.B.Eds. Springer US,49–68,1989.

[44]Beccari,C.V.;Casciola,G.;Romani,L.Non-uniform interpolatory curve subdivision with edge parameters built upon compactly supported fundamental splines. BIT Numerical Mathematics Vol.51,No.4,781–808, 2011.

[45]Conti,C.;Gemignani,L.;Romani,L.Exponential pseudo splines:Looking beyond exponential B splines. Journal of Mathematical Analysis and Applications Vol.439,No.1,32–56,2016.

[46]Botsch, M;Sorkine, O.On linear variational surface deformation methods.IEEE Transactions on Visualization and Computer Graphics Vol.14,No.1, 213–230,2008.

[47]Terzopoulos,D.;Platt,J.;Barr,A.;Fleischer, K.Elastically deformable models.ACM SIGGRAPH Computer Graphics Vol.21,No.4,205–214,1987.

[48]Davies,R.H.;Cootes,T.F.;Taylor,C.J.A minimum description length approach to statistical shape modeling.In:Information Processing in Medical Imaging.Insana,M.F.;Leahy,R.M.Eds.Springer Berlin Heidelberg,50–63,2001.

[49]Styner,M.A.;Rajamani,K.T.;Nolte,L.P.;Zsemlye, G.;Sz′ekely,G.;Taylor,C.J.;Davies,R.H.Evaluation of 3D correspondence methods for model building.In: Information Processing in Medical Imaging.Taylor, C.;Noble,J.A.Eds.Springer Berlin Heidelberg,63–75,2003.

[50]Aguet,F.;Jacob,M.;Unser,M.Three dimensional feature detection using optimal steerable filters.In: Proceedings of the IEEE International Conference on Image Processing,II-1158–II-1161,2005.

[51]Badoual,A.;Schmitter,D.;Unser,M.An inner product calculus for periodic functions and curves. IEEE Signal Processing Letters Vol.23,No.6,878–882,2016.

[52]Jacob,M.;Blu,T.;Unser,M.An exact method for computing the area moments of wavelet and spline curves.IEEE Transactions on Pattern Analysis and Machine Intelligence Vol.23,No.6,633–642,2001.

P.Garc′?a-Amorena was ranked in the top ten of the Spanish regional university entrance exam of Catalunya in 2009 and obtained double-bachelor degrees in mathematics and industrial engineering at Barcelona Tech(UPC), Spain.He received an Introduction to Research Grant from the Spanish National Research Council(CSIC)to work on 3D image reconstruction from deformation scans.In 2014 he won a La Caixa award for top undergraduate Spanish students. In 2016 he obtained his master degree in computational science and engineering at EPFL,Switzerland.His research interests include numerical methods and mathematical foundations for shape modeling,computer vision,and image processing. and Instrumentation Program,National Institutes of Health,Bethesda,USA,conducting research on bioimaging. Dr. Unser has held the position of associate Editor-in-Chief(2003–2005)for the IEEE Transactions on Medical Imaging.He is currently a member of the editorial boards of SIAM J.Imaging Sciences,IEEE J.Selected Topics in Signal Processing,and Foundations and Trends in Signal Processing. He is the founding chair of the technical committee on Bio Imaging and Signal Processing(BISP) of the IEEE Signal Processing Society.Prof.Unser is a Fellow of the IEEE(1999),an EURASIP Fellow(2009),and a member of the Swiss Academy of Engineering Sciences. He is the recipient of several international prizes including three IEEE-SPS Best Paper Awards and two Technical Achievement Awards from the IEEE(2008 SPS and EMBS 2010).

M.Unser is the professor and director of EPFL’s Biomedical Imaging Group, Lausanne,Switzerland.His primary area of investigation is biomedical image processing.He is internationally recognized for his research contributions to sampling theory,wavelets,the use of splines for image processing,stochastic processes,and computationalbioimaging.He has published over 250 journal papers on those topics. He is the author with P.Tafti of the book An Introduction to Sparse Stochastic Processes,Cambridge University Press,2014. From 1985 to 1997,he was with the Biomedical Engineering

Open Access The articles published in this journal are distributed under the terms of the Creative Commons Attribution 4.0 International License(http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use,distribution,and reproduction in any medium,provided you give appropriate credit to the original author(s)and the source,provide a link to the Creative Commons license,and indicate if changes were made.

Other papers from this open access journalare available free of charge from http://www.springer.com/journal/41095. To submit a manuscript,please go to https://www. editorialmanager.com/cvmj.

his master degree in bioengineering and biomedical technologies from the′Ecole Polytechnique F′ed′erale de Lausanne (EPFL),Switzerland, in 2013.He was with the Advanced Clinical Imaging Technology Group,Siemens, at the Center for Biomedical Imaging, Switzerland,where he was one of the main contributors working on brain-imaging software and related imageprocessing algorithms.Currently,he is a Ph.D.student at the Biomedical Imaging Group,EPFL,where he is working on spline-based shape representation and segmentation problems. He has developed several segmentation and tracking methods in the field of biomedical imaging.

1 Biomedical Imaging Group, ′Ecole Polytechnique F′ed′erale de Lausanne (EPFL), 1015 Lausanne, Switzerland. E-mail: D. Schmitter, daniel. schmitter@epfl.ch();P.Garc′?a-Amorena,pablo. garcia-amorenagarcia@epfl.ch; M. Unser, michael. unser@epfl.ch.

Manuscript

2017-01-24;accepted:2017-04-23