Zhen Chen

Adobe, USA Danny M Kaufman
email: dannykaufman@gmail.com

Mélina Skouras
email: melina.skouras@inria.fr

Etienne Vouga
email: etiennevouga@cs.utexas.edu

Complex Wrinkle Field Evolution

Keywords: Singularity Singularity Upsampling Upsampling Consistency Consistency Computing methodologies → Mesh geometry models discrete differential geometry, complex representation, wrinkle fields, singularities singularity Captures singularity

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Fig. 1. We propose Complex Wrinkle Fields (CWF s), a new discrete wrinkle model that enables the resolution of highly detailed wrinkle patterns on coarse base-mesh geometry. The CWF representation consists of a positive number a per vertex encoding the wrinkle amplitude, a one-form ω per edge to model wrinkle frequency, and a complex number z per vertex to represent wrinkle phase, coupled via a weak variational consistency condition ensuring that z can capture singularities while also being as compatible with ω as possible ( § 3.1). We equip the CWF representation with a novel temporal interpolation algorithm ( § 4) and a spatial upsampling method ( § 5) that together allow for smooth interpolation between wrinkle patterns represented on surfaces by CWF s (leftmost and rightmost column), and base-mesh-independent rendering of arbitrarily high-resolution wrinkle patterns. Together these contributions make it possible to smoothly evolve wrinkle patterns between two prescribed keyframes (middle columns) with automatic merging, splitting, and reconnection of wrinkles as necessary via smooth sliding of singularities across the surface (zoomed-in figures in middle columns). Please check the Interpolation Results extra video (00:38-00:53; note this is an additional video separate from the main supplemental video) for the corresponding wrinkle animation.

We propose a new approach for representing wrinkles, designed to capture complex and detailed wrinkle behavior on coarse triangle meshes, called Complex Wrinkle Fields. Complex Wrinkle Fields consist of an almosteverywhere-unit complex-valued phase function over the surface; a frequency one-form; and an amplitude scalar, with a soft compatibility condition coupling the frequency and phase. We develop algorithms for interpolating between two such wrinkle fields, for visualizing them as displacements of a Loop-subdivided refinement of the base mesh, and for making smooth local edits to the wrinkle amplitude, frequency, and/or orientation. These algorithms make it possible, for the first time, to create and edit animations of wrinkles on triangle meshes that are smooth in space, evolve smoothly through time, include singularities along with their complex interactions, and that represent frequencies far finer than the surface resolution.

INTRODUCTION

Across widely ranging spatial and temporal scales, wrinkles on surfaces are a fundamental geometric structure. We encounter these structures in our daily interactions with the moving folds and creases in cloth, skin and films, and likewise, we regularly observe them in natural phenomena such as the slow evolution of sand dunes and the rapid rippling of shallow water. Wrinkles on surfaces thus critically enrich otherwise coarse geometric structures with important visual details, while manufactured wrinkle patterns enable the design of complex structures and material behaviors [Chen et al. 2021;[START_REF] Lähner | DeepWrinkles: Accurate and Realistic Clothing Modeling[END_REF][START_REF] Zuenko | Wrinkles, Folds, Creases, Buckles: Small-Scale Surface Deformations as Periodic Functions on 3D Meshes[END_REF]. The real-world examples illustrated in Figure 2 demonstrate the characteristic and highly complex features of surface wrinkling behaviors in the wild. Locally, wrinkles are generally characterized by a dominant frequency and amplitude, both of which vary inhomogeneously but smoothly over the surface. At the same time singularities punctuate wrinkle patterns (see e.g. circled regions in Figure 2). These singularities are the branch points of the phase field where multiple peaks and troughs join together. Near these singularities wrinkles are fundamentally non-bandlimited: as we approach a singularity, circulating an arbitrarily small distance around it results in an arbitrarily large shift in phase. Figure 2, bottom row, shows several frames of wrinkle evolution during twisting of a sleeve of a cloth. Singularities appear, merge, and slide over the surface as wrinkles shift, split, and reconnect during their evolution.

In the smooth setting (see § 2.1) wrinkles have an appealing and compact representation as amplitude and phase signals-onsurface; on the other hand the above characteristic features pose significant modeling challenges when it comes to discretizing and evolving wrinkles on triangle meshes. In this work we address these challenges.

Contributions

First, because the evolving singularities and high-frequency wrinkles in fine wrinkle patterns are generally not captured at practical mesh resolutions, we construct a discretization of wrinkles called the Complex Wrinkle Field (CWF) that captures in-element singularities and multiple sub-element wrinkling periods ( § 3). The CWF representation enables mesh-independent resolution of fine wrinkles.

Second, to enable smooth design, editing, control, animation, and evolution between wrinkle patterns we construct a mechanics-based algorithm for continuous and temporally-coherent interpolation between arbitrary wrinkle patterns on surfaces ( § 4). Our algorithm takes as input two CWF endpoints ("keyframes") and solves a boundary value problem approximating a geodesic path in the space of inextensible shells. The results are complex, shifting wrinkle patterns, where singularities evolve automatically to support the necessary branching and merging of wrinkles needed to interpolate between the keyframed patterns.

Third, to complement CWF's ability to model sub-element wrinkle resolution, we derive a new subdivision method that maps CWFs on triangle meshes to amplitude and phase on vertices generated by Loop iterates ( § 5). This method resolves arbitrarily fine wrinkled geometry, including singularities, while avoiding artifacts seen in existing baselines (see § 5.3) and so provides high-quality upsampling of wrinkles for rendering and other downstream applications.

We extensively evaluate each of these contributions relative to state-of-the-art alternatives, and then demonstrate their joint application to animating wrinkle evolution on triangle surfaces, and to smooth, user-in-the-loop wrinkle editing (including both local and global manipulations). Across these applications we show that our contributions account for evolving singularities while generating smooth and intuitive interpolants between wrinkle patterns ( § 6). For example, in Figure 3 (bottom row), we show that, with proposed techniques, we can generate wrinkles with the double frequency and half of amplitude of the simulated results from [Chen et al. 2021]-effectively replacing thick cloth with thinner, silkier material without re-simulation. Furthermore, we can swipe through a smooth interpolation of the wrinkle patterns to select the desired look.

Despite the prominence of evolving wrinkles in the natural world, we emphasize that to date there is no model or method that supports temporal wrinkle interpolation on surfaces (in the sense of computing a path of spatially-and temporally-continuous wrinkle patterns between prescribed keyframes). Our CWF representation and algorithms are the first designed from the ground up to address the challenges imposed by evolving singularities during wrinkle interpolation. These challenges are generally unavoidable, even when the underlying surface and the prescribed keyframes are simple. For example, consider the torus interpolation in Figure 3 (top row) and in the supplementary Interpolation Results video at 00:02-00:19. Although neither keyframe contains singularities, it is topologically impossible to follow a smooth path between them without introducing (and then annihilating) singularities along the way.

In summary, Complex Wrinkle Fields (CWF) make it possible, for the first time, across all wrinkle inputs, to create, edit, and animate detailed wrinkles on triangle meshes that are smooth in space and evolve smoothly through time; while, at the same time, implicitly including the complex interactions of singularities, and resolving wrinkle frequencies far finer than the surface mesh resolution. Raw data and a reference implementation of our algorithms, as of the time of publication, can be found in the ACM Digital Library. The latest version of the code can also be found on GitHub1 .

PRELIMINARIES AND RELATED WORK 2.1 Wrinkles on Surfaces

Smooth Setting. A useful [START_REF] Knöppel | Globally Optimal Direction Fields[END_REF][START_REF] Knöppel | Stripe Patterns on Surfaces[END_REF] representation of wrinkles on smooth manifolds M is by a single complex smooth function z : M → C, where the wrinkle amplitude a is encoded as |z|, and the associated phase as θ = arg(z). Here we assume wrinkles are normal displacements of M with magnitude a cos θ = ℜ(z). The smoothness of z guarantees a smooth nonnegative amplitude, and a smooth associated phase, except at a set of singular points where z = 0 and along branch cuts between them. At these singular points the phase is undefined and the amplitude vanishes. Away from these singularities, the derivative of phase yields the (one-form-valued) wave frequency ω = dθ . It is important to note that the frequency field singularities are not the usual ones studied in geometry processing where ω = 0; rather, here ∥ω ∥ → ∞ as one approaches a singularity. Discrete Setting. Frequency and amplitude are the semanticallymeaningful parameters for encoding the wrinkling phenomena discussed in §1 and shown in Figure 2. Our goal then is to translate the above spectral representation of wrinkles from the smooth setting to triangle meshes T = {V , E, F }. As we motivated in the last section, a practical discretization of z on a triangulation by a set of discrete parameters u should satisfy three basic desiderata:

(1) The parameters u must allow for singularities; i.e., it should be possible to reproduce wrinkle patterns with branch points as illustrated in Figure 2. ( 2) High frequency wrinkling should be representable independent of the resolution of the underlying triangle mesh, and in particular, wrinkles should be allowed much more than a single period within a triangle. This decoupling (the main motivation for a spectral representation of wrinkles) is needed both because fine wrinkle patterns are common in the real world (as in the skin and cloth images in Figure 2) and because high frequencies are unavoidable in the vicinity of singularities.

(3) We must be able to render a smoothly wrinkled, high-resolution surface for all valid u, and this surface geometry must change continuously for continuous changes in u. This requirement represents a minimum foundation for doing temporally-coherent wrinkle interpolation.

Surprisingly, given how straightforward it is to represent wrinkles in the smooth setting, the most direct strategies for discretizing z (e.g., storing a complex number z i on each vertex) do not meet the above requirements, as we discuss shortly; constructing the right discrete representation to support temporal interpolation is a challenge in its own right. After reviewing related work we will present our discrete CWF model in §3, and compare with some natural alternative choices of discretization and demonstrate their shortcomings.

Prior Work

Spectral Representations of Wrinkles on Surfaces. Prior work has recognized the significant challenges in representing wrinkles on triangle meshes, and has offered a range of partial solutions. A common approach [START_REF] Aharoni | The smectic order of wrinkles[END_REF]Chen et al. 2021;[START_REF] Paulsen | Curvature-induced stiffness and the spatial variation of wavelength in wrinkled sheets[END_REF]] is to design for the wave frequency without regard to phase, using vector-field design tools from geometry processing; when phase is required (e.g., for rendering), it is recovered as a post-process by solving for the phase whose gradient most closely matches the frequency: either via least-squares optimization, or using techniques from the parameterization literature [START_REF] Diamanti | Integrable PolyVector Fields[END_REF][START_REF] Ling | Spectral Quadrangulation with Feature Curve Alignment and Element Size Control[END_REF][START_REF] Ray | Periodic Global Parameterization[END_REF][START_REF] Zhang | A Wave-Based Anisotropic Quadrangulation Method[END_REF]]. [START_REF] Knöppel | Stripe Patterns on Surfaces[END_REF] propose a particularly powerful variation of this idea: their algorithm solves for a stripe pattern (phase field) over arbitrary surfaces given a user-provided vector field, where singularities are placed at the centers of the triangles as needed to produce topological dislocations in the pattern. [START_REF] Noma | Fast Editing of Singularities in Field-Aligned Stripe Patterns[END_REF] build on this technique by allowing users to interactively edit the location of singularities in stripe patterns. Beyond designing stripe patterns, Knöppel et al.'s algorithm has been applied for real-time surface parameterization [START_REF] Lichtenberg | Real-time field aligned stripe patterns[END_REF]] and continuous fiber design [START_REF] Boddeti | Optimal design and manufacture of variable stiffness laminated continuous fiber reinforced composites[END_REF].

While these "frequency-first" strategies can produce beautiful, high-frequency wrinkle patterns on surfaces, they all suffer from a key limitation: smooth change in frequency does not yield a smooth change in phase. These strategies are therefore unsuitable for smooth animation of wrinkles due to temporal incoherence: they will produce a temporally discontinuous phase sequence, even when the input frequency sequence undergoes a smooth temporal evolution. For example, in the main supplementary video (02:59-03:44) and the Comparisons video (00:05-00:59), we show an animation of rotating waves on a torus, where we apply the method described by Chen et al. [2021] and [START_REF] Knöppel | Stripe Patterns on Surfaces[END_REF] to every frame. Although the individual frames are spatially smooth, they are not temporally coherent.

Vector Field Design. Over decades, many methods have been developed for vector field design. [START_REF] Vaxman | Directional Field Synthesis, Design, and Processing[END_REF] and de Goes et al. [2016a] summarize the classical ways to design and discretize vector fields on triangle meshes. In general, there are several ways to represent discrete vector fields. One popular approach specifies each vector with respect to a local Cartesian or polar coordinate system [START_REF] Diamanti | Designing N-PolyVector Fields with Complex Polynomials[END_REF][START_REF] Knöppel | Globally Optimal Direction Fields[END_REF]]. This representation is agnostic to singularities and well-suited for methods that automatically place them [START_REF] Bommes | Mixed-Integer Quadrangulation[END_REF][START_REF] Diamanti | Integrable PolyVector Fields[END_REF][START_REF] Wenzel | Instant Field-Aligned Meshes[END_REF][START_REF] Panozzo | Frame Fields: Anisotropic and Non-Orthogonal Cross Fields[END_REF][START_REF] Ray | Geometry-Aware Direction Field Processing[END_REF]]. However, enforcing global integrability of the vector field in this representation is difficult and requires non-linear constraints or integer variables. Another approach represents vector fields implicitly in terms rotation angles on dual edges [START_REF] Crane | Trivial Connections on Discrete Surfaces[END_REF][START_REF] Li | Representing Higher-Order Singularities in Vector Fields on Piecewise Linear Surfaces[END_REF][START_REF] Ray | N-Symmetry Direction Field Design[END_REF]]. This representation is particularly useful for parameterizing vector fields while maintaining full control over singularity placement and index [START_REF] Fisher | Design of Tangent Vector Fields[END_REF][START_REF] Solomon | Optimal Transport-Based Polar Interpolation of Directional Fields[END_REF][START_REF] Zhang | Vector Field Design on Surfaces[END_REF]] but integrability is even more elusive: two levels of integrability constraints are required (one for integrating the angles to well-defined vectors and another for integrating the vectors to a well-defined scalar field). As we discuss in § 3, neither of these representations are suited for representing temporally-coherent, high-frequency wrinkle fields with moving singularities on discrete triangle meshes.

Vector Field Subdivision. Subdivision operators have been designed that extend Loop subdivision to one-forms on triangle mesh edges [de Goes et al. 2016b;[START_REF] Wang | Edge Subdivision Schemes and the Construction of Smooth Vector Fields[END_REF]]. We will use this scheme as a part of our approach ( § 5.1). [START_REF] Custers | Subdivision Directional Fields[END_REF] designed a subdivsion scheme for tangent directional fields on triangle faces. However, these methods do not cover the corresponding subdivision of amplitude and phase, which are essential for upsampling wrinkles.

Vector Field Temporal Interpolation. Temporally interpolating vector fields has long been studied in geometry processing. [START_REF] Zavala-Hidalgo | A new temporal interpolation method for high-frequency vector wind fields[END_REF] interpolate high-frequency vector wind fields by decomposing these fields into a basis of empirically-derived orthogonal functions, and then interpolating eigenmodes in time. [START_REF] Chen | Design of 2D Time-Varying Vector Fields[END_REF] propose a framework for designing time-varying vector fields by solving a spatial-temporal Poisson problem; their method implicitly generates bifurcations and singularities in the vector field. [START_REF] Sato | Editing Fluid Animation Using Flow Interpolation[END_REF] propose an algorithm that takes two entire time sequences of velocity fields and produces time sequences that blend between them. [START_REF] Solomon | Optimal Transport-Based Polar Interpolation of Directional Fields[END_REF] use optimal transport theory to match the singularities in two input vector fields, which allows for their direct interpolation. However, vector field interpolation, on its own, does not solve the problem of interpolating wrinkles, where the phase and amplitude must also be interpolated, and, as discussed above, there is no obvious way to recover temporally-coherent phase from frequency. Spectral Wrinkle Evolution. There are several methods for simulating spectrally-represented wrinkles forwards through time from an initial state: in a line of work, Jeschke, Wojtan, and collaborators [2018;2015;2017] propose several approaches to efficiently evolve water waves based on Airy wave theory [1842], which uses a sum of sinusoidal functions to linearly approximate the motion of surface waves on the body of water. [START_REF] Zuenko | Wrinkles, Folds, Creases, Buckles: Small-Scale Surface Deformations as Periodic Functions on 3D Meshes[END_REF] simulate wrinkle motion on 3D surfaces via reaction-diffusion [START_REF] Turk | Generating Textures on Arbitrary Surfaces Using Reaction-Diffusion[END_REF][START_REF] Witkin | Reaction-Diffusion Textures[END_REF]. Although these approaches generate animations of wrinkles given initial conditions (initial value problem), they do not address the wrinkle keyframe interpolation problem, which is completely different (boundary value problem).

Shell Geodesics. [START_REF] Heeren | Time-Discrete Geodesics in the Space of Shells[END_REF] propose a way to find the geodesic between two thin shells based on the minimization of strain dissipation. This idea is further explored for time discretization of geodesic calculus on certain Riemannian manifolds [START_REF] Rumpf | Variational time discretization of geodesic calculus[END_REF], understanding the geometry of the space of shells [START_REF] Heeren | Exploring the Geometry of the Space of Shells[END_REF], or finding geodesics between two general elastic shapes [START_REF] Ezuz | Elastic Correspondence between Triangle Meshes[END_REF][START_REF] Sassen | Nonlinear Deformation Synthesis via Sparse Principal Geodesic Analysis[END_REF]. [START_REF] Von-Tycowicz | Real-Time Nonlinear Shape Interpolation[END_REF] adopt a similar idea and design a special subspace optimization problem to interpolate between elastic body poses in real time. Our wrinkle interpolation algorithm follows in the footsteps of these methods.

Surface Augmentation. Methods have also been proposed to augment coarse surfaces with fine wrinkles as an alternative to expensive, fine-scale physical simulation. Indirect techniques have been proposed for upsampling cloth [START_REF] Bergou | TRACKS: Toward Directable Thin Shells[END_REF][START_REF] Wang | GPU-Based Simulation of Cloth Wrinkles at Submillimeter Levels[END_REF] and skin [START_REF] Rémillard | Embedded Thin Shells for Wrinkle Simulation[END_REF] by simulating a higher-resolution mesh that is constrained to stay close to an existing, coarse simulation. [START_REF] Rohmer | Animation Wrinkling: Augmenting Coarse Cloth Simulations with Realistic-Looking Wrinkles[END_REF] and [START_REF] Gillette | Real-Time Dynamic Wrinkling of Coarse Animated Cloth[END_REF] add fine wrinkle detail to coarse cloth simulations by tracing compression direction fields in a temporal coherent manner. [START_REF] Cutler | An Art-Directed Wrinkle System for CG Character Clothing[END_REF] allow users to design wrinkle patterns on a set of reference poses, then add weighted combinations of the designed wrinkles to each frame of simulation based on local stresses. Data-driven methods for enriching coarse meshes with fine detail, for example, by representing wrinkles as normal maps [START_REF] Lähner | DeepWrinkles: Accurate and Realistic Clothing Modeling[END_REF] or displacements [START_REF] Chen | Synthesizing cloth wrinkles by CNN-based geometry image superresolution[END_REF][Chen et al. , 2021;;[START_REF] Santesteban | LearningâĂŘBased Animation of Clothing for Virtual TryâĂŘOn[END_REF][START_REF] Zhang | Deep detail enhancement for any garment[END_REF], are likewise an active and rapidly evolving area of research but require significant preprocessing and generally depend on a corpus of data.

COMPLEX WRINKLE FIELD DISCRETIZATION

In the smooth setting, a single complex function z suffices to represent a wrinkle pattern. Unfortunately, z admits no straightforward discretization that safisfies the properties listed in § 2.1. In this section, we introduce our complex wrinkle field (CWF) representation, the design decisions that motivate it, and some of the pitfalls with alternative discretizations. In the following sections we will then discuss how to both (a) temporally interpolate between CWFs and (b) upsample (and so render) CWFs while preserving these properties.

High-level Summary. There are three key design decisions that motivate CWFs:

• To represent high-frequency wrinkles on a mesh, unlike in the smooth setting, a phase-based representation is insufficient.

The discretization must additionally track frequency. • For the wrinkle representation to capture singularities, the frequency field cannot be globally integrable in the usual sense of being the derivative of a real-valued function, and cannot even be locally integrable near singularities. We therefore eschew integrability entirely and allow arbitrary frequencies. • Since frequency is not integrable, we cannot and should not exactly enforce that frequency is the derivative of phase. But we also do not want phase and frequency to be incompatible. We therefore require soft compatibility of phase and frequency, in a way that is well-defined at singularities and enforces d arg(z) ≈ ω away from singularities.

A CWF therefore consists of three sets of degrees of freedom:

(1) a one-form ω jk on the edges of T (representing the wrinkle pattern frequency); ( 2) a real scalar a j per mesh vertex (representing the wrinkle amplitude); and (3) a complex number zj per vertex (encoding the wrinkle phase). There are no constraints on the ω jk and a j , but z is required to be approximately unit and to satisfy a weak compatibility condition with respect to ω. The reader immediately interested in the formal definition of a CWF may skip to § 3.1; in what follows we motivate and justify each of the above decisions.

Need for Both Frequency and Phase. The most direct way to discretize z is as a discrete function on the mesh vertices. The strength of this approach is that it uniquely pins down the wrinkle amplitude |z| and wrinkle phase arg(z) at each vertex. However, z alone is unable to represent wrinkles with high frequencies, as illustrated in Figure 4: at most one wrinkle period per edge is possible, since the phase change arg(z j ) -arg(z i ) across the edge is bounded by 2π . By contrast, a frequency one-form ω can encode arbitrarily high frequencies independent of mesh resolution; however, a discrete representation based on frequency alone is unsuitable for applications like ours requiring temporal coherence, since wrinkle phase is underdetermined given frequency alone.

A mixed representation that tracks both phase and frequency simultaneously combines their advantages, and eliminates their disadvantages. This observation motivates a discrete representation that includes both a real-valued one-form ω i j on the mesh edges, and a complex number z j per mesh vertex. (In the CWF representation, we further decompose z j into parts as z j = a j zj with a j , zj as our DOFs rather than z j itself, for reasons we cover below.)

No Integrability Constraints on Frequency. In the smooth setting, frequency is defined as the derivative of phase. A natural question when discretizing frequency is whether to require it, by analogy, to be (discretely) integrable. Several important subtleties make doing so impractical. Note that, in the smooth setting, ω is not globally integrable in the usual sense: ω is the exterior derivative of a section of an S 1 -bundle over M-with-singularities-removed. The line integral of ω around a singularity does not necessarily vanish and is instead a multiple of 2π . Requiring ω = 0 around every closed curve in M would constrain the wrinkle pattern to have no singularities. Similarly, in the discrete setting, if a triangle jkℓ contains a singularity, the discrete curl on that triangle ω jk + ω k ℓ + ω ℓj does not vanish. Rather than selectively enforcing integrability of ω on only some triangles (which would require explicitly tracking which triangles contain singularities of which index using integer variables, etc.) we instead follow in the footsteps of previous work [START_REF] Knöppel | Stripe Patterns on Surfaces[END_REF][START_REF] Zuenko | Wrinkles, Folds, Creases, Buckles: Small-Scale Surface Deformations as Periodic Functions on 3D Meshes[END_REF] and allow arbitrary ω. The singularities (d): the result when we apply our subdivision algorithm from § 5 to visualize a CWF with z j = exp(3π ix j ) and the unrelated, incompatible

ω jk = (0, 3π ) • (v k -v j ).
Although the rendered wrinkled mesh is smooth, it is riddled with mesh-dependent singularities, and ω has little semantic relationship to the waves in the resulting pattern. We insist on soft compatibility between ω and z to avoid such artifacts. (e): this time, we use

ω jk = (3π, 0) • (v k -v j ),
a compatible frequency that gives us a valid CWF . When rendered, this discrete CWF is indistinguishable from the exact pattern.

in the wrinkle pattern are then implicitly determined by where ω fails to be integrable.

Frequency-Phase Consistency. Since we do not require ω to be integrable, ω cannot in general be exactly the derivative of phase. Recall, however, that the motivation for including frequency in our discrete representation is so that it can supply the information missing from z j (alone) about the jump in phase across a mesh edge. For this purpose, frequency is only useful if, far away from singularities, we do have ω ≈ d arg(z). (See Figure 4d for an example where frequency and phase are totally uncorrelated.) We therefore require the phase variables in a CWF to be "as compatible as possible" with the frequency, (a) given that frequency is not necessarily integrable and (b) accounting for the necessary presence of singularities.

In the smooth setting, one formulation of the above soft compatibility condition is to require z to be compatible in a least-squares sense by minimizing an energy similar to M d arg(z) -ω 2 dA.

(1)

We cannot use this energy as-is, as it is formally undefined near singularities and at branch cuts. To avoid the discontinuity in arg(z) at branch cuts, we rewrite the constraint in terms of z = z/|z| [START_REF] Knöppel | Stripe Patterns on Surfaces[END_REF],

z = arg min z M ∥(d -iω) z∥ 2 dA (2) s.t . | z| = 1,
which has the added benefit of eliminating some of the nullspace from Equation (1) (which is invariant under rescalings of z).

Singularity Handling. Equation ( 2) remains ill-posed at singularities (where neither z nor ω are well-defined) and ill-conditioned close to the singularities, where ∥ω ∥ → ∞.

There are many potential strategies for relaxing Equation ( 2) to regularize the behavior at singularities. Our approach is to relax the unit-norm constraint, so as to allow z to vanish at singularities, using a Ginzburg-Landau-type term [START_REF] Robert | Energy-driven pattern formation[END_REF][START_REF] Viertel | An Approach to Quad Meshing Based on Harmonic Cross-Valued Maps and the Ginzburg-Landau Theory[END_REF]:

z ∈ Opt ω = arg min z M ∥(d -iω) z∥ 2 + z 2 -1 2 dA. (3)
For later reference, we define:

E compat ( z, ω) = M ∥(d -iω) z∥ 2 dA, E unit ( z) = M z 2 -1 2 dA. ( 4
)
We next discretize this variational characterization of z to define the final, discrete soft compatibility conditions we impose to our CWF representation.

Complex Wrinkle Field Representation

Motivated by the above discussion, we define CWF to consist of the following degrees of freedom:

(1) a (not necessarily integrable) frequency one-form, encoded as a real number ω jk on each directed edge of the mesh (with ω jk = -ω k j ); (2) a real number a j and complex number zj on each vertex of the base triangle mesh T , which together encode the approximate amplitude and direction of z j = a j zj discretizing the smooth setting's amplitude-phase field z;

(3) with the z j constrained to be softly compatible with the frequency one-form, d arg(z) ≈ ω, via discretization of Equation (3) on the base triangle mesh.

To discretize the first, least-squares, compatability term of Equation (3), we apply the formulation proposed by [START_REF] Knöppel | Stripe Patterns on Surfaces[END_REF], as it is designed to be robust in the case that ω jk is larger than 2π on an edge jk,

Êcompat ( z, ω) = edges k ℓ A k ℓ zk exp(iω k ℓ ) -zℓ 2 , (5)
and directly apply piecewise-constant discretization for the second, unit-length penalty term,

Êunit ( z) = vertices j A j zj 2 -1

INTERPOLATION

In this section we present an algorithm for temporal interpolation of CWFs: given two CWFs as boundary conditions, ω 0 , a 0 , z0 and

ω 1 , a 1 , z1 , we seek a smooth interpolant γ (t ) = (ω[t], a[t], z[t]
) through the space of CWFs which matches the prescribed boundary conditions at t = 0, 1 and which approximates the behavior of real-world wrinkled materials deforming over time.

The challenge with computing such an interpolation is that while there are many smooth choices of interpolant γ (t ), they are often qualitatively unacceptable. For physically-based wrinkle evolution, wrinkles should slide over the surface rather than disappear in one location and then reappear at a target location, and singularities should slide smoothly over the surface. Otherwise plausible-seeming interpolants that are not specifically designed to promote physical fairness of the wrinkle motion regularly produce obvious visual artifacts that are unacceptable for many applications such as animation or editing. For many examples please see § 4.3 and the Interpolation Results video.

Our Approach

To compute physics-inspired, smooth paths between boundary conditions in the space of CWFs, we begin with the closely-related problem of computing geodesics in the space of thin shells (see e.g. the work of [START_REF] Heeren | Time-Discrete Geodesics in the Space of Shells[END_REF]). For our application, the key idea from this work boils down to defining a metric on shell space that measures the length of a path between configurations in terms of the integrated norm of strain-rate along that path. A shortest geodesic in this space is then, effectively, the path that requires doing the least work on the system when deforming the shell along that path.

To compute comparable wrinkling geodesics for CWFs we measure a bending strain 2 given in terms of the combined geometry of the underlying wrinkle-free base surface M and the wrinkle displacements [Chen et al. 2021],

ϵ = I -1 II -Ī I -a cos θω T ω (8) = I -1 II -Ī I -ℜ(z)[d arg z] T [d arg z] , (9)
where I and II are the first and second fundamental forms of the base surface and Ī I is the second fundamental form of the shell's assumed rest state. Following [START_REF] Heeren | Time-Discrete Geodesics in the Space of Shells[END_REF], we can then define our distance on paths γ (t ) by

d (γ ) = 1 0 M ∥ ε[z(t )]∥ 2 dA dt, (10)
where, assuming a hyperelastic homogeneous and isotropic material, we apply the StVK material norm for ∥ • ∥; please see the supplemental document ( § 1) for details.

Computing Geodesics. Geodesics are then constrained minimizers γ of Equation ( 10) satisfying z(t ) ∈ Opt ω (t ) . However, computing these geodesics by directly minimizing Equation ( 10) is difficult. In the supplemental document ( § 1.2) we derive, after a sequence of approximations and additional simplifying assumptions, an approximate solution to Equation (10) given by

a(t ) = (1 -t )a 0 + ta 1 (11) z(t ) = arg min E smooth ( z) + cE opt ( z, ω) (12)
and an explicit formula for ω (t ) given in Equation ( 34) ( § 1.2 in the supplemental document), where

E smooth ( z) = 1 2 M д bd 1 0 ˙z 2 dtdA (13) д bd = a 0 2 ω 0 4 + a 1 2 ω 1 4 2 ( 14
)
approximates Equation (10) after a and ω are substituted in and simplified-effectively giving a weighted Dirichlet energy in spacetime promoting smoothness of wrinkle phase; and

E opt ( z, ω) = 1 0 E compat ( z[t], ω[t]) + E unit ( z[t]) dt (15)
provides a penalty term, with stiffness c, for the constraint z ∈ Opt ω that z and ω are approximately compatible. We use c = 10 3 д ave in all of our examples, where

д ave = 1 0 M a 2 ∥ω ∥ 4 dA dt ( 16
)
is the spacetime averaged amplitude squared times frequency quartic, in order to make the cE opt unit consistent. If c is too low, the interpolant γ (t ) fails to satisfy the CWF compatibility constraint, yielding noisy wrinkle pattern evolution and many extra singularities. If c is too high, the optimization problem (12) becomes numerically unstable and does not converge to a smooth solution.

Please see § 7.2 in the supplemental document for some experiments probing the effect of the choice of c on our results and the suitability of our empirically selected value.

Discretization. Given two CWF boundary conditions ω 0 , a 0 , z0 and ω 1 , a 1 , z1 and a desired number of intermediate frames N , we discretize time into N steps of size δt = 1/N and minimize a discretization of Equation ( 12) over the base triangle mesh T with a globalized Newton solve. Please see § 1.3 in the supplemental document for additional details, including the full formulas for all our discrete analogues of the above energy terms used in our final solve. Within our Newton implementation we use SuiteSparse's parallel implementation of supernodal sparse Cholesky decomposition [START_REF] Chen | Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate[END_REF]] as our linear solver, and terminate when converged to an objective gradient norm smaller than 10 -6 .

Linear Substepping. The cost of computing N interpolating frames using our algorithm is dominated by the cost of solving Equation ( 12), which scales roughly linearly in N . When many frames are needed, for example when rendering videos, this computation is wasteful, since consecutive frames will be almost identical. For large N , we suggest using our algorithm to compute N ′ = N /k guide frames instead, and then linearly interpolating the CWF variables for k steps between each guide frame. As a rule of thumb, linear interpolation is an adequate approximation to a geodesic path between guide frames if the distance traveled by wrinkles during that time does not exceed the average mesh edge length. We use N = 200, with N ′ = 50 and k = 4 in all of our videos, and provide timing information for the examples in the paper in Table 1 (supplemental document). We further provide some experiments showing the effect of the choice of N ′ on the interpolant quality in the supplemental document ( § 7.1).

Handling Incompatible Boundary Conditions. The above algorithm assumes that provided boundary conditions are valid CWFs. If applied to boundary conditions that severely violate CWF soft compatibility, our interpolant will move abruptly near t = 0 and 1 (since Equation (15) will enforce the soft compatibility condition everywhere except at the infeasible, prescribed endpoints). Since in applications these boundary conditions will often be provided by the user (as keyframes for animation, e.g.) we propose pre-and post-processing steps that allows use of our algorithm even when the boundary conditions have incompatible frequency and phase.

In particular, notice that the discrete compatibility term in Equation (5) implies the constraint per edge zk exp(iω k ℓ ) = zℓ .

(

) 17
This allows us to decompose any frequency one-form ω into a consistent, constraint-satisfying portion, ω com , and a residual δω; that is, there is a unique one-form δω, with -π ≤ (δω) jk < π on each edge jk, such that ω com = ωδω and z exactly satisfy Equation ( 17) on every edge. To handle non-CWF boundary conditions, we decompose the frequencies ω 0 and ω 1 into their consistent and inconsistent parts, and apply interpolation to the CWFs boundary conditions (ω 0 com , a 0 , z0 ) and (ω 1 com , a 1 , z1 ) instead. We then interpolate δω using the same explicit interpolation applied for ω (t ) given in Equation (34) ( § 1.2 in the supplemental document), and add this frequency component back into the final interpolant:

γ (t ) = (ω (t ) + δω (t ), a(t ), z(t )).

Results

We demonstrate that our interpolation algorithm yields high-quality wrinkle evolution given a variety of keyframes pairs. Even more interpolation examples, with keyframes created by local edits to wrinkle patterns, are discussed in § 6. Important note: it is difficult to assess the quality of the interpolant, and especially temporal coherence, from a sparse set of stills. Please see the Comparisons supplemental video for animations of these examples.

Change in Wrinkle Frequency. Increasing the frequency of a wrinkle pattern on a compact surface like the fertility model (Figures 5 and6), though conceptually simple, involves a complex path through the space of wrinkle fields since it is impossible to have "fractional" numbers of wrinkles. New wrinkles must be born as the frequency increases, with singularities appearing in pairs and sliding over the surface to "unzip" a new wrinkle.

Change in Wrinkle Direction. We also show several examples of interpolation between wrinkle patterns that have been globally or locally rotated by ninety degrees relative to each other in Figures 1 and7. The amplitude and phase plots show singularities appearing, splitting, and merging in a complex dance in both cases.

Comparison to Other Approaches

We compare our method against some baselines.

Linear Interpolation. It is natural to try to temporally interpolate between CWF boundary conditions by linearly interpolating the constituent variables:

a(t ) = (1 -t )a 0 +ta 1 ; ω (t ) = (1 -t )ω 0 +tω 1 ; z(t ) = (1 -t ) z0 +t z1 .
Such an interpolant lacks physical meaning and moreover does not stay within the space of CWFs, since there is no reason to expect ω (t ) and z(t ) to stay approximately compatible during the interpolation. See Figure 10 and Figures 9, 10 in the supplemental document ( § 6) for some examples of artifacts that arise when using this linear scheme. Note that replacing linear interpolation of z and ω with a more complicated non-linear interpolant would not help resolve this incompatibility issue.

For comparison on a more didactic example, consider the simple case of a plane wave rotating by ninety degrees, z 0 = exp(ix ) and z 1 = exp(iy). On a sufficiently small patch centered at the origin, we expect interpolation to rotate the plane wave, with an interpolant resembling z(t ) = exp(i[x cos π t 2 + y sin π t 2 ]). In Figure 8 we see that when using our proposed interpolation algorithm, the interpolant is indeed rigid rotation of the wrinkles. On the other hand, linear interpolation shears the wrinkles, which midway through the interpolation degenerate along a singular line. [START_REF] Knöppel | Stripe Patterns on Surfaces[END_REF]. In their paper, Knöppel et al. propose a method to extract vertex phase information from a provided frequency one-form by solving an eigenvalue problem. They also propose a spatial interpolation scheme for extending phase into triangles by placing singularities at their centers as needed. While well-suited for generating individual, static geometries, Knöppel et al. 's method is not designed to give temporally coherent motion when run on a sequence of frames with smoothly-changing frequency. Indeed, if we apply their method to ω (t ) (computed using Equation (34) in the supplemental document § 1.2) to generate a phase interpolant z(t ), we get beautiful individual frames, as expected (Figure 10) but see that the frames are discontinuous in time (see the main supplemental video at timestamp 03:44-04:04).

Moreover, Knöppel et al. focus exclusively on phase fields and do not discuss wrinkles with amplitude. We make a good-faith attempt to incorporate amplitude in the method (See Figure 9) by linearly interpolating the keyframe vertex amplitudes in time and then barycentrically blending them into each triangle. The resulting artifacts are due to a mismatch between where Knöppel et al. places phase singularities (at centers of triangles) and where amplitude vanishes (which can only be at vertices).

Chen et al. [2021].

In their recent work, Chen et al. describe a method for rendering high-resolution wrinkled surfaces given frequency and amplitude fields on the surface. To compare against this work, we again use Equations ( 12), Equation (34) (in the supplemental document § 1.2) to compute a(t ) and ω (t ) and apply Chen et al. to each frame of the interpolation. See Figure 10 and the Comparisons supplemental video for the results. Unlike Knöppel et al.

[2015], Chen et al. 's method does consider wrinkle amplitude, and so can produce smooth wrinkled surface geometry. However, their approach again suffers from temporal incoherence. Moreover, we observe their method also can fail to capture the sliding of singularities over the surface, instead producing degenerate, noisy phase fields; see e.g., Figure 10 around t = 0.5.

Keyframes from Physical Simulation. Finally, we compare the behavior of our interpolant to the results of physically simulated wrinkle evolution. We start by simulating the twisting of a 1-meter-high (diameter 0.25 meters, thickness 0.318 mm, Young's Modulus 0.821 MPa, density 472.6 kg/m 3 , and Poisson's ratio 0.243), 88k-vertex cylinder mesh, by 80 degrees with C-IPC's [START_REF] Li | Codimensional Incremental Potential Contact[END_REF] opensource implementation of shell finite elements. The cylinder's top and bottom boundaries are pinned to rotate (time-stepped at 0.04 seconds/frame) in opposite directions at 10 degrees/s without gravity or external forces. To avoid simulating the transient dynamics before wrinkles appear, we initialize the cylinder slightly twisted (by 8 degrees in both direction) and run the simulation for 80 frames (so that cylinder ends are twisted by 40 degrees in both directions at the last frame). During this motion, wrinkles rotate and merge; see Figure 11, top, and the Comparisons supplemental video.

To compare these simulated results with corresponding intermediate frames computed using CWF interpolation, we take the first and last frame of the simulation as interpolation inputs. Specifically, we decimate the cylinder mesh to 688 vertices and then apply a tension-field-theory-based static solver [START_REF] Skouras | Designing Inflatable Structures[END_REF]], implemented by Chen at al. [2021], to compute a wrinkle-free base mesh for every frame. We then manually set the frequency at the two endpoint keyframes by counting the number of wrinkles in the corresponding frames of the C-IPC simulation, and scale the wrinkle amplitude to preserve surface area of the cylinder. In Figure 11, bottom, we see that the CWF interpolated frames then well-align with the input keyframes and, in-between, obtain qualitatively similar merging and deformation to the simulation results. That said, these interpolated frames do have expected and notable differences from the simulated results: deformation paths differ (e.g., with wrinkles merging at different times) and, of course, the CWF model (as it considers only a single primary frequency) does not capture the secondary wrinkles observable in the simulation results.

WRINKLE UPSAMPLING

Representing fine wrinkle patterns on coarse meshes is only useful if there is a way to actually see them. Given a triangle mesh T and a CWF on that mesh, rendering the wrinkled surface requires both upsampling T itself (yielding a smooth canvas on which to place wrinkles, free of meshing artifacts) and turning the CWF into a displacement map that can be applied to the upsampled mesh.

To that end, we propose the following upsampling strategy:

Linearly blended amplitude Knöppel's phase Wrinkled mesh CWF wrinkled mesh Fig. 9. [START_REF] Knöppel | Stripe Patterns on Surfaces[END_REF]'s algorithm is designed only for stripe patterns (phase fields), not wrinkles with amplitude. Naively combining the phase field produced by Knöppel et al's method with linearly-interpolated amplitude (through barycentric blending from triangle corners) yields artifacts near wrinkle singularities (the red zoomed-in region). In this case, the coarse mesh has a constant amplitude on each vertex. Moreover, without our CWF subdivision algorithm, meshing artifacts are apparent in the rendered result (the black zoomed-in region).

• For upsampling T , we focus exclusively on Loop subdivision [1987], as it is the most commonly-used triangle mesh subdivision scheme. • In the remainder of this section, we will describe a rule for mapping per-vertex complex numbers z and per-edge frequencies ω from T k , the input mesh after k levels of Loop refinement, to T k +1 . The quantities on T 0 can be read directly from the CWF (z i = a i zi on the vertices), and then upsampled to any desired level of refinement. • To materialize the wrinkled mesh at subdivision level k, we apply the displacement |z i | ni cos arg z i to vertex V i , where ni is the vertex normal.

We use this approach to render all smooth wrinkled surfaces in the paper (we apply 4 or 5 levels of subdivision, depending on the frequency of wrinkles in each example). We have not analyzed this process to prove that the limit surfaces belong to some smoothness class, though in practice all of our results are visually smooth, including near singularities. Although the above procedure is designed to render CWFs, it can also be used to visualize z, ω pairs that do not obey approximate compatibility. We have observed that the resulting upsampled wrinkles are also smooth, although they include extraneous singularities not reflected in z or ω (see Figure 4d).

Preliminaries. We now describe the upsampling maps z k → z k +1 and ω k → ω k +1 . To avoid cumbersome superscripts we will write T for T k and T ′ for T k+1 , and similarly use primes to denote quantities on T k +1 . We write V for the vertices of T and S 0 for the Loop subdivision mask, so that V ′ j = [S 0 V ] j .

Upsampling ω

Research in geometry processing [de Goes et al. 2016b;[START_REF] Wang | Edge Subdivision Schemes and the Construction of Smooth Vector Fields[END_REF]] has established a subdivision mask S 1 on one-forms that is compatible with the Loop subdivision mask and the discrete differential operator, in the sense that for any function f on T , S 1 df = dS 0 f . On closed meshes, we use this operator directly, and set ω ′ = S 1 ω.

On meshes with boundary, we modified the standard Loop formula and de Goes et al. 's formula slightly to keep boundary corners sharp (see § 2 of the supplemental document for details).

Upsampling z

There are two obvious ideas for how to upsample the z values: either by subdividing z's real and imaginary parts separately, or its magnitude and phase. Neither approach works, for essentially the same reasons discussed in § 3 as motivating the use of frequencies in addition to phase in the CWF representation: in regions where the frequency is high, phase values at vertices undersample the wrinkles. Computing z ′ by averaging these aliased samples yields an equally-aliased z ′ . We show a comparison between our approach (described below) and these two naive baselines in Figure 12.

Our Approach.

To upsample z without aliasing the wrinkles, we need to incorporate frequency information from ω into the upsampling procedure. We build up to the final formulas for computing z ′ on the vertices of T ′ by first examining some simpler special cases. We will use the following assumption as our guiding principle for deriving our subdivision rules: in a local patch of the surface, the phase of z might change quickly, but the frequency does not. This assumption allows to extrapolate and interpolate frequency values away from the edges on which they are defined.

Upsampling on a Line Segment. Let us first consider the case of a line segment with vertices V 0 and V 1 , whose z-values are z 0 and z 1 , and whose edge frequency is ω. Suppose we need to estimate z at a point P = (1-α )V 0 +αV 1 . Since we assume frequency changes slowly, we can extrapolate z at P from z 0 by z(P ) ≈ z 0 exp(iαω); similarly from z 1 by z 1 exp(i (1-α )ω). These formula do not necessarily agree unless |z 0 | = |z 1 | and ω and z are exactly compatible. We reconcile the two votes by barycentrically blending them:

z(P ) = (1 -α )z 0 exp(iαω) + αz 1 exp(i (1 -α )ω). (18)
Notice that z(P ) interpolates z 0 and z 1 .

Upsampling in a Triangle. We use a similar strategy for evaluating z at any point P = 3 j=0 α j V j within a triangle {V 0 , V 1 , V 2 }, whose corresponding z values are z 0 , z 1 and z 2 , and whose edges have frequencies ω 01 , ω 12 , and ω 20 . Similar to the line segment case, we can estimate z(P ) by extrapolating from a triangle corner V j by z(P ) ≈ z j exp iω j→P , where

ω j→P := α j+1 ω j (j+1) + α j+2 ω j (j+2)
is the frequency on the segment V j P, as measured at vertex V j using ω j (j+1) and ω j (j+2) as a basis for evaluating ω in any direction. Barycentrically blending this formula from three corners gives us:

z(P ) = 3 j=1 α j z j exp i α j+1 ω j (j+1) + α j+2 ω j (j+2) . (19)
This triangle interpolant is similar to the one proposed by [START_REF] Zuenko | Wrinkles, Folds, Creases, Buckles: Small-Scale Surface Deformations as Periodic Functions on 3D Meshes[END_REF] with one key difference: they formulate Equation (19) purely for phase and propose barycentrically interpolating amplitude separately. The Zuenko and Harders approach is unsuitable for triangles containing singularities, since upsampled amplitude does not neccessarily vanish at the singularites, creating visual artifacts (Figure 16).

Loop Subdivision Rules for z. We now generalize Equation (19) to compute z ′

i on the Loop-upsampled vertex V ′ i of T ′ . We must consider two cases: V ′ i might be a repositioned vertex already present in T (an even vertex), or an entirely new vertex that came from splitting an edge of T (an odd vertex). We separately describe how to compute z ′ in each case; see Figure 13 for a visual summary of all the rules. Odd Rule. For an odd vertex, the standard Loop mask is

V ′ j = 1 8 (V 0 + V 1 ) + 3 8 (V 2 + V 3 )
on a "diamond" with common edge V 2 V 3 . This mask can be rewitten

V ′ j = 1 2 P 0 + 1 2 P 1 P 0 = 1 4 V 0 + 3 8 V 2 + 3 8 V 3 P 1 = 1 4 V 1 + 3 8 V 2 + 3 8 V 3 . (20)
Notice that these formula combine barycentric sampling within a triangle (to compute P 0 and P 1 ) followed by sampling the midpoint of a line segment (to compute V ′ j .) We can convert this interpolation into two single triangle interpolations at P 0 , w.r.t. triangle {V 0 , V 2 , V 3 } and P 1 , w.r.t. triangle {V 1 , V 2 , V 3 }, together with line segment interpolation at V ′ j , w.r.t. {P 0 , P 1 }. We can transform these rules for computing V ′ j into rules for computing z ′ j by applying the formulas from the special cases discussed above: we use Equation ( 19) to compute z at P 0 and P 1 . Next, we would like to use Equation ( 18) to compute z ′ j . There are two obstacles: first, we need the frequency evaluated on the line segment P 0 P 1 . Second, this line segment generally does not lie on T , but rather cuts through ambient space, whereas frequencies are intrinsic to T 's cotangent space.

We therefore modify the approach of Equation ( 18) to account for these differences. First, we sample ω at P 0 by interpolating the values of ω at the edges of triangle V 0 V 2 V 3 into the interior (see § 3 of the supplemental document for more details; in particular note that ω (P 0 ) is a cotangent vector and not a scalar). We use this frequency to extrapolate z(P 0 ) to V ′ j along the intrinsic edge from P 0 to P 1 :

z P 0 →V ′ j = z(P 0 ) exp iω (P 0 )(P * 1 -P 0 )/2 , (21)
where P * 1 is the location of P 1 after rigidly unfolding triangle V 1 V 2 V 3 along the common edge V 2 V 3 to lie in triangle V 0 V 2 V 3 's tangent plane. We compute z P 1 →V ′ j analogously. Then the final z ′ j is the average of the corresponding contributions from P 0 and P 1 :

z ′ j = 1 2 (z P 0 →V ′ j + z P 1 →V ′ j ). (22)
Even Rule. Next, we consider an even vertex V ′ 0 corresponding to V 0 on T , with neighbors V 1 , . . . , V n . V ′ 0 has position

V ′ 0 = (1 -αn)V 0 + n j=1 αV j , ( 23
)
where α is the traditional even-neighbor Loop weight and depends on the valence n of V 0 [START_REF] Loop | Smooth Subdivision Surfaces Based on Triangles[END_REF]]. As in the odd case, we can refactor Equation ( 23) as a convex combination of points linearly interpolated within a single triangle:

V ′ 0 = n-1 j=0 1 n P j , P j = (1 -2γ )V 0 + γV j+1 + γV j+2 , (24)
for γ = nα 2 . We now follow the same recipe as in the odd case: we compute z(P j ) using Equation ( 19), and write z ′ 0 as the average of contributions extrapolated from each P j . 18) in our upsampling algorithm with some alternatives, for a CWF on the Stanford bunny. From left to right: the nine-parameter Clough Tocher cubic scheme [START_REF] Farin | Triangular Bernstein-BÃľzier patches[END_REF]], the side-vertex scheme used by [START_REF] Jeschke | Water Wave Animation via Wavefront Parameter Interpolation[END_REF], and Equation (19). Only Equation (19) successfully resolves the phase behavior near singularities (zoomed-in regions).

Other Upsampling Methods. In Figure 16, we compare our upsampling algorithm against two baselines from the literature [Chen et al. 2021;[START_REF] Zuenko | Wrinkles, Folds, Creases, Buckles: Small-Scale Surface Deformations as Periodic Functions on 3D Meshes[END_REF]. Zuenko and Harders render fine wrinkles on coarse meshes by first linearly subdividing the mesh, interpolating phase onto the fine mesh using a variant of Equation ( 19) and amplitude using barycentric interpolation, and then smoothing the result with SPHERIGON [START_REF] Volino | The SPHERIGON: a simple polygon patch for smoothing quickly your polygonal meshes[END_REF]] as a post-process to remove some of the artifacts from the wrinkled mesh geometry. Since they separately upsample wrinkle amplitude and phase, there is no guarantee that amplitude vanishes at phase singularities, leading to noticable noise (in the red and black zoomed-in regions, for instance). Moreover the use of linear subdivision on T means that the edges of the coarse mesh remain noticeable in the final rendered image, despite the post-processing.

Like in our approach, Chen et al. [2021] use Loop subdivision to upsample T , as well as the wrinkle amplitude and phase, but their method requires that the input frequency and phase are exactly compatible. This restriction makes their approach unsuitable for rendering CWFs. In Figure 16 we make a best-effort comparison by projecting ω to the closest globally-integrable frequency using mixed-integer programming, but doing so severely distorts the wrinkle pattern, and many singularities are missing from the final rendered image.

WRINKLE DESIGN

In this section, we demonstrate applications of our CWF algorithms for user-based design and editing of wrinkled surfaces via interpolation. We first develop tools for the creation and editing of CWF keyframes and then show their application. Please also see our main and Interpolation Results videos for more details and results.

Adding Wrinkles to Surfaces

For adding wrinkles to a surface a user first provides a mesh T = {V , E, F } and selects k vertices as source points s (green points in Figure 17a) and at each of those points a desired wrinkle frequency vector v s (the yellow arrow). A user then also specifies a region of influence in s ⊂ V for each source point (the white and red regions [2021]'s upsampling scheme only works for curl-free frequency fields. To use their method, we project ω onto the space of globally integrable oneforms, but doing so distorts the wrinkle pattern, which is missing many singularities.

Designed sources

Extended (a, ω ) Wrinkled surface Amplitude Phase Fig. 17. Adding wrinkles to the Spot cow. From left to right: The user specifies desired wrinkle directions (yellow vectors) and regions of influence for each chosen direction (red regions). The user also assigns an amplitude to each region. We extend a and ω from the user input to the entire mesh. Finally, we solve for z and apply our upsampling algorithm to materialize high-resolution wrinkled geometry. We show the amplitude and phase of the final wrinkle pattern, after upsampling, on the right.

in Figure 17a)3 and a target amplitude a s for the wrinkles within the region of influence. Let S be the set of all source points {s i } k i=1 . We call the set of all vertices in any region of influence the influenced vertices in V = in s i . The influenced faces in F are those with at least one influenced vertex, and likewise for the influenced edges in E .

From this input we solve for a wrinkle pattern on T . The process can be divided into two parts: (1) input extension, where we extend the user-provided amplitude and frequency samples to an a and ω on the entire mesh, and ( 2) wrinkle synthesis, where we compute z from a and ω and use the upsampling technique introduced in § 5 to export the final wrinkled geometry.

We first apply the vector heat method [START_REF] Sharp | The Vector Heat Method[END_REF] to extend the user-provided frequencies v s i to a vertex-based vector field v i on V . Then for any influenced edge E i j ∈ in E we set

ω i j = 1 2 (v i • E i j + v j • E i j ),
and for any edge that is not influenced, we set ω i j = 0. We compute per-vertex amplitudes by solving

arg min a E (a) s.t.      a i = 0, V i in V a s i = c s i , ∀s i ∈ S (25)
where

E (a) = 1 2 a T La + k i=1 V j ∈in s i (a j -a s i ) 2 A j + F ℓmn ∈in F a 2 ℓ + a 2 m + a 2 n 3 (ω ℓm + ω mn + ω nℓ ) 2
and L is the positive semi-definite cotangent Laplacian matrix, A j is the vertex barycentric area, and ω ℓm + ω mn + ω nℓ is the discrete curl. The first term promotes smoothness of a; the second term measures agreement with the user-specified amplitudes in the affected regions, and the final term couples the a to the ω to ensure that amplitude goes to 0 at the singularities of ω. Finally, for the wrinkle synthesis step, we solve the eigenvalue problem proposed by [START_REF] Knöppel | Stripe Patterns on Surfaces[END_REF] to determine z and then directly apply the subdivision algorithm from § 5 to (ω, a, z).

Remarks. The procedure above is used to create local wrinkle patterns. If global patterns are desired instead (such as the wrinkle patterns on the bunny in Figure 7), we first compute ω using [START_REF] Knöppel | Globally Optimal Direction Fields[END_REF]'s algorithm, and globally scale the frequency and set a global constant amplitude to achieve the desired wrinkle characteristics. Then we find the consistent z values on the vertices and upsample the result using the wrinkle field synthesis algorithm described above.

Editing Wrinkles

Next, we describe a tool for editing wrinkles, given an existing mesh T and CWF ( z0 , a 0 , ω 0 ). There are four steps in the editing pipeline:

• The user picks a region of the mesh to edit (red in Figure 7).

• The user edits the wrinkles in that region. Options we implemented are: rotating the frequency, changing the frequency magnitude, and changing the wrinkle amplitude. • We blend the new wrinkle pattern in the edited region with the existing one inside an interface region around the chosen region (green region in Figure 7). Wrinkles away from the edited region and its interface do not change. • We interpolate between the original CWF and the edited one, using our interpolation algorithm.

The third step requires some elaboration. Let F edit ⊆ F be the set of faces where the user has chosen to make edits, and let ω 1 and a 1 be the requested frequency and amplitude in that region. We need to construct a new z1 , which satisfies (1) z1 = z0 far away from F edit ; and (2) z1 is consistent with ω 1 within F edit . To avoid a discontinuity in frequency and amplitude on the edited surface, we create an interface region F inter of user-specified size around F edit . The remaining mesh faces will not change; we call them F fixed .

Let E edit be the edges of the triangles in F edit , and E fixed and V fixed the edges and vertices of F fixed , respectively. The above requirements can be written as follows:

• z1 j = z0 j , for any V j ∈ V fixed . • | z1 j | = 1, for any V j ∈ V \ V fixed . • z1 k exp(iω k j ) -z1 j = 0 for any edge E jk ∈ E edit .
• The interface frequency and amplitude fields are smooth.

We first interpolate ω 1 into the interface region. To get a frequency field that is as smooth as possible, we minimize the Dirichlet energy of that field with Dirichlet boundary conditions:

arg min ω F ℓmn (ω ℓm + ω mn + ω nℓ ) 2 + V m E mn c mn ω mn 2 s.t        ω jk = ω 0 jk , E jk ∈ E fixed ω jk = ω 1 jk , E jk ∈ E edit , (26)
where c mn is the cotangent weight w.r.t. edge E mn , and the first and second term are the discrete curl and divergence respectively. Once we have the frequency field, the amplitude field is the optimal solution of the following modified Dirichlet energy:

arg min a 1 2 a T La + F ℓmn a 2 ℓ + a 2 m + a 2 n 3 (ω ℓm + ω mn + ω nℓ ) 2 s.t      a j = a 0 j , V j ∈ V fixed a j = a 1 j , V j ∈ V edit , (27)
which asks for an amplitude that is as smooth as possible while ensuring (via the first term) that the amplitude vanishes at singularities.

Here L is the positive semi-definite cotangent Laplacian. Finally, we compute an approximately-consistent z inside the interface and edited regions by solving an optimization problem generalizing the strategy from [START_REF] Knöppel | Stripe Patterns on Surfaces[END_REF]:

min ũ,s E compat (P ũ + su 0 , ω) s.t. P ũ + su 0 2 M 0 = 1 (28)
where ũ is a vector of size |V | -|V fixed | representing the unknown values of z in the non-fixed region, u 0 j = z0 j for any V j ∈ V fixed (with other entries zero), M 0 is the mesh barycentric mass matrix, and P is the inclusion matrix that extends ũ to the whole domain by inserting zeros for vertices in the fixed region. The single scalar value s allows for rescaling (but no other changes) within the fixed region. The inclusion of s in the optimization is required since the more straightforward restriction of [START_REF] Knöppel | Stripe Patterns on Surfaces[END_REF]'s method to the non-fixed region,

min ũ E compat (P ũ + u 0 , ω) s.t. P ũ + u 0 2 M 0 = 1, (29)
enforces an average scaling of ũ within the non-fixed region separate from, and generally inconsistent with, the scaling of u 0 in the fixed region.To transform Equation ( 28) so that it is once again an eigenvalue problem, we make the substitutions v = [ ũ, s] and M ′ 0 = P T M 0 P. From ũ we compute u = P ũ/s + u 0 and then set z1 j = u j /|u j |. Notice that (ω 1 , a 1 , z1 ) is not necessarily a valid CWF since z1 may not satisfy the soft compatibility condition exactly (and it is not generally possible to satisfy them without changing the zj within the fixed region). We can render the wrinkle field nonetheless, and use it as a keyframe for wrinkle evolution by projecting it to a CWF using the pre-processing described in § 4.

Results

We perform a variety of local edits to wrinkles on surfaces and visualize the wrinkle evolution during the editing process. Figures 6 and 19 illustrate local frequency increase, with a corresponding decrease in amplitude in the former example. In Figure 7 a user demands a local rotation of the wrinkles by 90 degrees, while maintaining the frequency magnitude and amplitude. We also support different types of edits on different patches. In Figure 18, we enlarge the frequency of the patch on Spot's body by 2.5 times, while at the same time rotating the direction of wrinkles on the head by 90 degrees. These examples show that we can interpolate between extreme frequency or direction differences, and the singularities nevertheless slide smoothly over the surface without temporal discontinuities or other artifacts seen in results from previous work.

Please see the Interpolation Results supplemental video, and the supplemental document ( § 5), for more examples of editing wrinkles, including some examples of modifying the output of wrinkles computed via physical simulation.

LIMITATIONS AND FUTURE WORK

In this paper we made several design decisions that could be fruitfully revisited in future work in order to generalize CWFs and our interpolation or upsampling algorithms, such as restricting wrinkles to waves with a single (but spatially-varying) frequency, rather than a superposition of frequencies [START_REF] Rémillard | Embedded Thin Shells for Wrinkle Simulation[END_REF]; and limiting refinement to Loop subdivision. The comparison to wrinkled materials (e.g., the simulated cloth in Figure 11) indicates significant potential value in extending this work to model additional, secondary, higher-frequency wrinkle fields superimposed upon the first. Another particularly interesting direction of future work is to study the use of CWFs for solving PDEs involving waves other than the interpolation boundary value problem: for example, we believe that it's likely that cloth dynamics would benefit from CWF kinematics as it would allow simulation of wrinkles that have much higher frequency than the base mesh. Finally, complex wrinkle fields could be used to model skin wrinkles or fingerprints; though as argued by [START_REF] Zuenko | Wrinkles, Folds, Creases, Buckles: Small-Scale Surface Deformations as Periodic Functions on 3D Meshes[END_REF], high-quality skin wrinkles would likely require extending CWFs to waves with non-sinusoidal profile.

Interpolation Performance. Interpolating CWFs is currently far from real-time (please refer the time table in the supplemental document for more details); our implementation took 4 minutes to compute keyframe interpolation on average, with 32 seconds minimum and 12 minutes maximum among all the examples. One problem is that our optimization problem is not a convex optimization: although the other two terms (E smooth and E compat ) in Equation ( 12) are quadratic, E unit is nonlinear and not convex, and this nonlinearity is enough to cause Newton's method difficulty. Performance could be improved by refactoring the optimization problem (including perhaps by improving the initial guess) or parallelizing the linear solver on the GPU. Note that the size of the linear system (in terms of number of non-zeroes) grows linearily in the number of requested frames N . To improve the efficiency, we proposed the use of guide frames ( § 4.1). We study the performance vs. quality tradeoff of the choice of number of guide frames in the supplemental document ( § 7.1); careful tuning of this parameter could achieve a bettter tradeoff than our (conservative) use of N ′ = 50.

Preventing Collisions. When amplitudes are large, it is possible for neighboring wrinkle periods to collide, or for wrinkles to "poke through" distant portions of the base mesh. Our work currently does not detect or prevent such collisions in any way.

Alternative Incompatibility Norm. When defining Opt ω , we choose to use the L2 norm to measure the incompatibllity between ω and z (Equation ( 2)). Given that singularies are usually sparsely distributed on a wrinkled surface, a promising alternative is the L1 norm, which promotes such sparsity, but is more difficult to optimize.

More General Wrinkle Parameterization. Several of our editing examples take as input wrinkle geometry computed by spectral wrinkle simulators [Chen et al. 2021;[START_REF] Zuenko | Wrinkles, Folds, Creases, Buckles: Small-Scale Surface Deformations as Periodic Functions on 3D Meshes[END_REF]. In these cases we could directly convert the output of the simulators into CWF variables. For more general wrinkled surfaces (computed by a finite element simulator [START_REF] Narain | Adaptive Anisotropic Remeshing for Cloth Simulation[END_REF]], for instance), a missing step is how to decompose a high-resolution, wrinkled mesh into the coarse but smooth base mesh T and the CWF variables. It's unclear how to best perform this decomposition, especially when wrinkles are complex, aliased, and contain many singularities. In addition to optimization-based approaches, another promising avenue is training a neural network to perform the decomposition. Improved Wrinkle Physics. Our interpolation is physics-inspired, but due to simplifying assumptions in the derivation such as neglecting the inertia of wrinkles and the effect of base mesh stretching on the wrinkle evolution, our interpolant is not guaranteed to exhibit fully realistic wrinkle dynamics. One possible improvement is to rederive the interpolation formula using the full reduced-order elastic energy derived by Chen et al. [2021] instead of assuming inextensible shells. Additionally, the current model interpolates the CWF variables completely independently of any movement or deformation of the base mesh. A scheme that jointly optimizes for a coupled base mesh and CWF interpolant could result in wrinkle evolution with higher physical realism. Along similar lines, our current model does not consider inhomogeneities or special features of the base mesh, such as hems or seams, which significantly affect fabric wrinkling in the real world. A possible direction for further exploration is to extend our work by imposing spatial boundary conditions on the CWF variables at the locations of stiff seams.

Fig. 2 .

 2 Photo courtesy Flickr user: nevil zaveri Photo courtesy Flickr user: Chris

Fig. 3 .

 3 Fig.3. Two examples of CWF interpolation between keyframes specified at t = 0 and 1. Top row: a torus, with wrinkles keyframed to rotate by ninety degrees; bottom row: a dress example fromChen et al. [2021], where we emulate replacing the cloth with a thinner material by uniformly doubling the wrinkle frequency and halving the amplitude. For the rotating torus example, although neither keyframe contains singularities, it is topologically impossible to follow a smooth path between them without introducing (and then annihilating) singularities along the way. Here our CWF interpolation allows wrinkles to break apart before reconnecting again. The dress demonstrates the possibility of using our algorithms to reasonably edit physical wrinkles and generate corresponding smooth animation. Videos of these two examples can be found in the Interpolation Results supplementary video at 00:02-00:19 (for the torus) and 03:47-04:03 (for the dress).

Fig. 4 .

 4 Fig. 4. We show different discretizations of a plane wave associated with z = exp(3π ix ). (a): the exact wrinkle pattern z , visualized as a normal displacement with magnitude ℜ(z ). (b): we sample z onto a coarse mesh and plot arg(z j ); the discrete phase is badly aliased. (c): attempting to visualize the wrinkles by displacing each vertex by ℜ(z j ) yields a surface that bears no resemblance to the expected wrinkle pattern in subfigure (a).(d): the result when we apply our subdivision algorithm from § 5 to visualize a CWF with z j = exp(3π ix j ) and the unrelated, incompatible ω jk = (0, 3π ) • (v k -v j ). Although the rendered wrinkled mesh is smooth, it is riddled with mesh-dependent singularities, and ω has little semantic relationship to the waves in the resulting pattern. We insist on soft compatibility between ω and z to avoid such artifacts. (e): this time, we use ω jk = (3π, 0) • (v k -v j ), a compatible frequency that gives us a valid CWF . When rendered, this discrete CWF is indistinguishable from the exact pattern.

Fig. 8 .

 8 Fig. 8. Failure of linear interpolation given boundary conditions z 0 = exp ix and z 1 = exp iy. We show the phase patterns at several intermediate frames using linear interpolation (top row) and our proposed interpolant (bottom row). Notice that linear interpolation introduces severe artifacts in the wrinkle pattern: wrinkles shear and degenerate along a sharp line midway through the interpolation.

Fig. 5 .Fig. 6 .

 56 Fig. 5. CWF interpolation between two wrinkle patterns on the fertility model, where the second keyframe (last column) has triple the frequency and one-third the amplitude of the first keyframe (first column). We show, from top to bottom, the amplitude, phase, and rendered wrinkles for several intermediate frames during interpolation. On the amplitude plot, blue indicates |z | = 0 (singularities) and red indicates larger amplitude. For the animation, please check the Interpolation Results video in the supplementary material at timestamp 01:46-02:01.

Fig. 7 .

 7 Fig. 7. CWF interpolation of two wrinkle patterns on the Stanford bunny, where the user has locally rotated the frequency field by ninety degree (red region of the top mesh). Top row: the user-specified keyframes. Second row: the upsampled amplitudes |z | = a | z | of the CWF , ranging from blue at singularities (|z | = 0) to red where the amplitude attains its maximum. Third row: the upsampled phase of the CWF . Bottom row: the rendered wrinkle patterns. We see nontrivial wave mergers via singularities sliding. Please check the Interpolation Results video in the supplementary material at timestamp 01:29-01:45 for the corresponding animation.

Fig. 10 .

 10 Fig. 10. We compare our interpolation algorithm to several baselines on an example of wrinkles rotating by ninety degrees on the Stanford bunny. Linear interpolation leaves the space of CWF s and consequently produces severe artifacts. Both Knöppel et al. [2015] and Chen et al. [2021] produce temporally incoherent results, though these discontinuities are hard to see in static frames and much more obvious in the main supplementary video (at 03:44-04:04) and the Comparisons supplementary video (at 01:00-01:56). Additionally, Chen et al. 's phase field becomes noisy midway through the interpolation. Note that the Chen et al. phase field differs from the others at t = 1 since their method recomputed phase from the provided frequency input.

Fig. 11 .

 11 Fig.11. We compare CWF interpolation to the results of physically simulated wrinkle evolution. We simulate the twisting of a cylinder by 64 degrees (from 16 to 80 degrees) with C-IPC[START_REF] Li | Codimensional Incremental Potential Contact[END_REF]] (top) and compare its results with CWF interpolation using the the first and last simulated frame (bottom) as boundary conditions. CWF smoothly interpolates the given keyframes with qualitatively-similar wrinkle merging and deformation (though it does not capture secondary wrinkles that emerge during simulation). See the Comparisons supplemental video (at 02:57-03:13) for an animated comparison.

  Fig. 12. A comparison of our z → z ′ upsampling algorithm (right) to naively applying the Loop subdivision mask to z 's amplitude and phase (left), and to z 's real and imaginary parts (middle). Neither baseline preserves the frequency of the wrinkles in the original CWF .

Fig. 15 .

 15 Fig. 15. A comparison where we replace Equation (18) in our upsampling algorithm with some alternatives, for a CWF on the Stanford bunny. From left to right: the nine-parameter Clough Tocher cubic scheme[START_REF] Farin | Triangular Bernstein-BÃľzier patches[END_REF]], the side-vertex scheme used by[START_REF] Jeschke | Water Wave Animation via Wavefront Parameter Interpolation[END_REF], and Equation (19). Only Equation (19) successfully resolves the phase behavior near singularities (zoomed-in regions).

Fig. 14 .

 14 Fig.14. We successively Loop-upsample the Stanford bunny, alongside the z and ω of a CWF on the bunny, and render the resulting displaced surface. After several rounds of subdivision, the wrinkle pattern appears and is stable under additional iterations of subdivision.

Fig. 18 .

 18 Fig. 18. CWF interpolation of two wrinkle patterns on Spot. Please refer to the Interpolation Results video in the supplementary materials for the corresponding wrinkle animation (timestamp 02:21-02:37).

Fig. 19 .

 19 Fig. 19. CWF interpolation of two wrinkle patterns on the Phantasma model, where the user has locally increased wrinkle frequency by 5×. Our algorithm successfully generates a smooth path between these challenging keyframes. For the corresponding wrinkle animation, please watch the Interpolation Results video at timestamp 00:54-01:11 in the supplementary materials.

ACM Trans. Graph., Vol. 1, No. 1, Article 1. Publication date: January 2023.

https://github.com/zhenchen-jay/Complex-Wrinkle-Field.git ACM Trans. Graph., Vol. 1, No. 1, Article 1. Publication date: January

,(6)where the A k ℓ and A j are edge and vertex barycentric area weights, respectively. We then requirez ∈ Ôpt ω = arg min z Êcompat ( z, ω) + Êunit ( z) ,(7)where Ôpt ω defines our discretized consistency conditions.Remarks. Note that in Equations (4) and (7) we have made an implicit choice to weight the compatibility and unit-norm terms equally. For regions away from the singularities, E compat ≈ 0, so that Equation (4) gives us unit z for any positive scaling of E unit . For regions near the singularities, E compat → ∞ unless z → 0 so that, for any scaling of E unit , the least-squares energy E compat is the dominant term, which ensures that z ≈ 0. These observations show that Opt ω is largely invariant to the relative scaling of the two constraint terms, and we make the simplest choice to weight them equally. Note also that Opt ω is a nontrivial subspace and contains more than a single element. At minimum, if z ∈ Opt ω , then so is its global phase shift exp(iθ ) z for every θ ∈ [0, 2π ).

We make the simplest assumption that the stretching strains in the wrinkled surface represented by a CWF are negligible. This decision is justified as compressive stresses should be mostly absorbed by wrinkling. That said, unless both CWF keyframes of an interpolation represent isometric deformations of the same underlying surface, there must be some amount of stretching involved and so it should be interesting to consider inclusion of stretching measures as well in future work. ACM Trans. Graph., Vol. 1, No. 1, Article 1. Publication date: January 2023.Complex Wrinkle Field Evolution • 1:7

ACM Trans. Graph., Vol. 1, No. 1, Article 1. Publication date: January 2023. Complex Wrinkle Field Evolution • 1:13

For simplicity we compute the region of influence as the n-ring neighborhood of the chosen points, for a user-specified n; one could also compute a geodesic disk.

ACKNOWLEDGMENTS

We thank Rahul Narain, Chris Wojtan, Justin Solomon, Max Wardetzky, Ulrich Pinkall, Keenan Crane, and Felix Knöppel for valuable discussions while preparing this work; Peter Schröder and Fernando de Goes for their help tracking down Loop one-form subdivision code; Minchen Li for guidance compiling the IPC codebase; and Kevin Song, Xinya Zhang, and Josh Vekhter for their support and advice. We thank the maintainers of the libigl, Polyscope, C-IPC, and Eigen libraries for contributing high-quality open-source software to our community, and the National Science Foundation (IIS-1910274) and Adobe Systems for their generous support.

Mesh Boundaries. For even vertices V 0 on the mesh boundary with neighbors V 1 and V 2 , we can once again write the Loop subdivision mask in terms of an average of linear interpolations along mesh edges, and apply the same recipe as in the even and odd interior cases above. In this case V ′ 0 = (P 0 + P 1 )/2, where

The odd boundary case involves sampling z on a mesh edge and so Equation ( 18) can be applied directly. V 0 ( 1 2 )

Fig. 13. Subdivision rules for even and odd interior and boundary vertices. First, sample ω and z (using Equation ( 19)) at the red points P i , using the barycentric weights specified in parentheses next to the vertices. Finally, extrapolate z (P i ) to the upsampled vertex (blue or yellow) using an analogous approach to Equation (21), and average these extrapolations using the weights in parentheses next to the P i .

Remark. Notice that the upsampling strategy explained above has the following useful property: it exactly reproduces plane waves, in the sense that if T is a piece of the plane, ω is the differential of a linear function θ (P ) = v • P, and z(= z) is consistent with ω, then under refinement the CWF (ω, 1, z) converges to (a phase shift of) the smooth plane wave ψ (P ) = cos(v • P ).

Results and Comparisons

We show successive upsampling of a CWF on the Stanford bunny in Figure 14. Wrinkles appear as soon as the subdivided mesh has high enough resolution to resolve their frequency. See all other figures of rendered wrinkles in the paper for more examples of our subdivision algorithm at work. Triangle Interpolation Alternatives. Several alternative formulas have been proposed for interpolating phase into a triangle from samples at its corners; in Figure 15 we compare our choice of Equation (19) to several potential alternatives: the side-vertex scheme used in Jeschke et al's work [2015], and the nine-parameter Clough-Tocher cubic scheme [START_REF] Farin | Triangular Bernstein-BÃľzier patches[END_REF]]. Notice that away from the singularities, all of these three approaches achieve reasonable results; however the side-vertex and Clough-Tocher cubic scheme fail to resolve neighborhoods of singularities. The reason for this failure is that both schemes express phase at interior points in terms of a rational or polynomial function of phase at the triangle vertices. It is mathematically impossible for such functions to produce singularities (where phase must have a branch point where it's undefined).