Prof. Loïc Barthe

Co-Supervisor

Prof Brian Wyvill

Prof Mathias Paulin

Marco Tarini

J P Lewis

Javier Von Der Pahlen

Miles Macklin

Sofien Bouaziz

Julien Valentin

Shahram Izadi

Daniel Rebain

Maurizio Kovacic

Wei Jiang

Kedi Xia

Sri Weiwei Sun

Raghu Malireddi

Abhishake Kumar Bojja

Nicolas Guillemot

Both man-made or natural objects contain repeated geometric elements that can be interpreted as primitive shapes. Plants, trees, living organisms or even crystals, showcase primitives that repeat themselves. Primitives are also commonly found in man-made environments because architects tend to reuse the same patterns over a building and typically employ simple shapes, such as rectangular windows and doors. During my PhD I studied geometric primitives from three points of view: their composition, simulation and autonomous discovery.

Résumé

Les objets naturels ou artificiels contiennent des éléments géométriques répétés pouvant être interprétés comme des formes primitives. Les plantes, les arbres, les organismes vivants ou même les cristaux représentent des primitives qui se répètent. Les primitives se retrouvent aussi couramment dans les environnements créés par l'homme, car les architectes ont tendance à réutiliser les mêmes modèles sur un bâtiment et utilisent généralement des formes simples, telles que des fenêtres et des portes rectangulaires. Au cours de ma thèse, j'ai étudié les primitives géométriques sous trois angles: leur composition, leur simulation et leur découverte autonome.

Dans la première partie, je présente une méthode pour retrouver la fonction par laquelle certaines primitives sont combinées. Notre système est basé sur un modèle de fonction de composition représenté par une surface paramétrique. La surface paramétrique est déformée via un alignement non rigide d'une surface qui, une fois convergée, représente l'opérateur souhaité. Cela permet de modéliser de manière interactive les opérateurs via une simple esquisse, ce qui résout un problème majeur de facilité de modélisation de la composition.

Dans la deuxième partie, je présente l'utilisation d'une nouvelle primitive pour les simulations de physique en temps réel. Cette primitive convient pour modéliser efficacement les déformations des cordes préservant le volume, mais également celles de structures plus complexes telles que les muscles. L'un des principaux avantages de notre approche est que notre primitive peut servir de représentation unifiée pour la détection de collision, la simulation et la déformation de peau.

Dans la troisième partie, je présente un cadre d'apprentissage en profondeur non supervisé pour apprendre et détecter les primitives. Dans un signal contenant une répétition d'éléments, le procédé est capable d'identifier automatiquement la structure de ces éléments (c'est-à-dire des primitives) avec une supervision minimale. Afin de former le réseau qui contient une opération non différenciable, un nouveau processus d'apprentissage en plusieurs étapes est présenté.

In the first part I present a method to reverse-engineer the function by which some primitives are combined. Our system is based on a composition function template that is represented by a parametric surface. The parametric surface is deformed via a non-rigid alignment of a surface that, once converged, represents the desired operator. This enables the interactive modeling of operators via a simple sketch, solving a major usability gap of composition modeling.

In the second part I introduce the use of a novel primitive for real-time physics simulations. This primitive is suitable to efficiently model volume-preserving deformations of rods but also of more complex structures such as muscles. One of the core advantages of our approach is that our primitive can serve as a unified representation to do collision detection, simulation, and surface skinning.

In the third part I present an unsupervised deep learning framework to learn and detect primitives.

In a signal containing a repetition of elements, the method is able to automatically identify the structure of these elements (i.e. primitives) with minimal supervision. In order to train the network that contains a non-differentiable operation, a novel multi-step training process is presented.

iii

List of Tables

Motivation

What is a geometric primitive? Looking up the word "primitive" in a dictionary does not give a clear idea. A primitive is defined as something "simple", "basic", "elemental", "unrefined", or "not derivative". But the definition that is the most relevant to my research can be found in [Collins, 2004]: "a form from which another is derived". Primitives are, in fact, defined not by their intrinsic properties but by their role within a larger structure. In this research, we consider geometric primitives to be recurring patterns or shapes that compose an object. Just looking through a window and observing how objects are composed should convince the reader that primitives are very common in our world. But the concept of geometric primitives is also important in computer graphics and computer vision, from geometric modeling to pattern recognition and others. We will then have a closer look at the role of primitives in these three domains: nature, computer graphics and computer vision.

Primitives in nature

Repetitions & tessellations. We find plenty of instances of repetitions of elements in nature.

Animals' skin often exhibit a recurring pattern. This is due to various reasons, involving chemistry [START_REF] Kondo | The reaction-diffusion system: a mechanism for autonomous pattern formation in the animal skin[END_REF] and natural selection [START_REF] Hugh | Adaptive coloration in animals[END_REF][START_REF] Endler | Natural selection on color patterns in poecilia reticulata[END_REF]. For example the skin of a snake is covered by small diamond-shaped scales; see Figure 1.1b. Similarly, many plants' structures are composed of the same primitives with very little variation stemming from the randomness of the natural process. For instance, pine cones are mainly composed of dozens of instances of a single primitive; see Figure 1.1a. The same thing can be observed for a spiral aloe or many other plants. The leaves of a tree can be seen as randomly varying instances of the same primitive.

Bees build hexagon-shaped structures called honeycombs in order to store larvae and honey; see Figure 1.1c. The striking coherence of this structure can be explained by the fact that this shape is a near-optimal solution to the problem of maximizing the cells volume for a given amount of wax [START_REF] Thomas | The honeycomb conjecture[END_REF][START_REF] Weaire | A counter-example to kelvin's conjecture on minimal surfaces[END_REF]. One can think of primitives as nature's response to an optimization problem. Because a shape is optimal for a certain function, it will naturally occur in many places.

Fractals. Many plants are not only structured repetitions of the same primitive but also selfsimilar at different scales. Shapes that are infinitely self-similar and have a fractal dimension are called fractals [START_REF] Benoit | The fractal geometry of nature[END_REF]. In nature however, infinite repetition is not attainable but there are many examples of few levels of recursion. For instance the Romanesco Broccoli can roughly be seen as a cone-shaped primitive; see Figure 1.2a. But this cone is filled with smaller cone-shaped elements. And these are themselves filled with similar elements. Trees often exhibit a fractal-like growth where the trunk splits into a few branches, each branch splits into smaller branches, etc; see Figure 1.2b. The fern leaves are composed of leaflets that are themselves composed of subleaflets and look similar to a whole leaf. All snowflakes are different but a lot of them exhibit some form of fractal-like shape. This is the case of the stellar dendrite in Figure 1.2c that is formed of three levels of hexagons.

Primitives in computer graphics

Efficient rendering. Three-dimensional objects are often represented by a mesh model that explicitly describes its surface [START_REF] Botsch | Projective dynamics: Fusing constraint projections for fast simulation[END_REF]. A mesh is a collection of vertices and polygonal faces. Faces are usually simple geometric primitives such as triangles or quadrilaterals. One of the reasons this representation is popular is that rendering triangles can be done very efficiently using a rasterizer [Akenine-Moller et al., 2018]. Rasterizers are typically implemented to run on GPUs that are designed to execute Single Instruction Multiple Data (SIMD). The OpenGL graphics API only offers to render a few primitives: points, lines, triangles, etc [START_REF] Shreiner | OpenGL programming guide: The Official guide to learning OpenGL[END_REF].

Figure 1.3: (left) A CSG composition tree. Primitives are combined pairwise using boolean operations. In the right branch, cylinders are combined by a union operator (∪ symbol). In the left branch, a cube and a sphere are combined by intersection (∩ symbol). At the root the final object is produced by computing the difference between two shapes (-symbol). (right) CAD softwares use CSG operations to model complex objects from simple primitives. (Sources: Zottie, Frank Hull) Composition modeling. In vector graphics, images are defined by combining together simple primitives (circles, polygons, curves, etc). These primitives can combined in many ways such as union, difference or intersection. This is exactly the same principle that is applied to 3D in constructive solid geometry (CSG) [START_REF] Requicha | Constructive Solid Geometry. TM[END_REF]; see Figure 1.3. Users can use simple primitives to model complex objects or scenes. L-systems use simple grammar rules to generate realistic shapes such as plants [START_REF] Prusinkiewicz | L-systems: from the theory to visual models of plants[END_REF][START_REF] Prusinkiewicz | The algorithmic beauty of plants[END_REF] or buildings [Parish and [START_REF] Yoav | Procedural modeling of cities[END_REF] with simple primitives.

Scene compression.

Streaming of 3D data is becoming popular and has applications in scientific visualization, cultural heritage, digital entertainment, healthcare, robotics, etc. One major challenge is the size of the data that needs to be transferred. Because bandwidth is limited it can cause latency issues which may not be acceptable for some applications, especially for ones where there is a real-time component involved. It is possible to use traditional compression methods but a more efficient way to compress scene data would be to leverage the repetitions and coherency in the geometry of our world. It is possible to represent scenes accurately with simple primitives e.g. man-made environments can often be described by a compact set of recurring elements, such as the windows of a building, or the arrangement of offices and desk in a university. These primitives can be detected or even be learnt and have a much smaller footprint than their corresponding point cloud [Kaiser et al., 2018;Tang et al., 2018a]. Physics simulation. Primitives are often used in physics simulation. In the finite element method (FEM) [START_REF] Cecil Zienkiewicz | The finite element method[END_REF], complex structures are decomposed into tetrahedra or hexahedra to homogenize the computations. In Lagrangian fluid simulation methods fluids are represented with particles that interact with each other and combine to represent a fluid [START_REF] Ihmsen | SPH Fluids in Computer Graphics[END_REF]. We can simulate complex objects by discretizing them with simple primitives. This is demonstrated by NVIDIA's FleX [NVIDIA, 2018], a popular physics simulation library based on position-based dynamics [START_REF] Müller | Position based dynamics[END_REF][START_REF] Macklin | Unified particle physics for real-time applications[END_REF]. In this framework, everything including soft bodies, cloth and fluids is represented as collections of particles connected by constraints. Using simple primitives also makes a physics engine simpler because there no need for specialized functions for each object. On GPU it is also more efficient to execute the same function on all the objects. Finally primitives are also used for collision detection [START_REF] Ericson | Real-time collision detection[END_REF].

Testing for collisions between two arbitrary shapes is expensive, so objects are often approximated by simpler primitives [START_REF] Ritter | An efficient bounding sphere[END_REF][START_REF] Larsen | Fast distance queries with rectangular swept sphere volumes[END_REF][START_REF] Rimon | Efficient distance computation using best ellipsoid fit[END_REF]. For instance the bounding boxes of two objects can be used [START_REF] Joseph | Finding minimal enclosing boxes[END_REF]. It is also common to use bounding volume hierarchies (BVH) [START_REF] Michael R Kaplan | The use of spatial coherence in ray tracing[END_REF][START_REF] Timothy | Ray tracing complex scenes[END_REF] where objects are roughly approximated by one primitive and if the collision test is positive the object's shape is refined into several primitives, which can themselves contain more levels of hierarchy.

Primitives in computer vision

Pattern recognition. The human visual system is built upon primitives. The cells connected to the retina respond to small circles within their respective regions. By measuring the response of cells of the visual cortex to visual stimuli, Hubel and Wiesel [1959] discovered that these cells build upon several dot detections to detect lines. They also identified complex cells that use groups of simple cells to detect the movement of lines [START_REF] David | Receptive fields, binocular interaction and functional architecture in the cat's visual cortex[END_REF]. Similarly, one could say that computer vision uses simple "primitives" to detect complex patterns and shapes [START_REF] Lecun | Deep learning[END_REF]. Nowadays convolutional neural networks (CNN) are the state of the art method for object detection tasks [START_REF] Liu | Deep learning for generic object detection: A survey[END_REF]. CNNs learn convolutional filters to detect patterns in images.

These filters detect recurring spatial patterns, i.e. primitives in the image domain, in pictures of objects they are trained to detect. On the first layer, the filters convolving on the input images are in the image domain. Typically the first layers' filters respond to low-level patterns such as edges and contrast. The ones in the hidden layers respond to primitives in the feature domain. The deeper they are, the higher-level patterns they will be looking for; see Figure 1.4. For instance for face detection, the next layer would be looking to find elements such as noses or eyes. While the last layer should be using looking for a arrangement of eyes, noses, mouth and ears to finally detect a face.

Scene simplification.

Primitives can also serve as geometric proxies to simplify a scene and make it easier to understand. In one of the earliest examples of this, Roberts [1963] attempted to reconstruct a three-dimensional object from a single photograph. His method was based on edge detection which were then converted to simple polygonal primitives to reconstruct the 3D object.

This method was only applicable to simple shapes with planar faces but more complex objects can be simplified as aggregates of simple primitives. Fischler and Elschlager [1973] proposed a system where objects like planes are abstracted by a template made of generalized cylinders. This template is not a fixed 3D model and has constrained parameters that vary to match images.

Similarly, in [Brooks et al., 1979] objects are decomposed as parts such as eyes and noses that are linked by spatial relations. This serves as a template that can be detected in images. It is a similar idea that is exploited in the recent work on "capsule networks" [START_REF] Sabour | Dynamic routing between capsules[END_REF][START_REF] Hinton | Matrix capsules with em routing[END_REF].

Body tracking. Body tracking has a wide variety of applications including computer-human interactions, performance capture, entertainment and telepresence. Thanks to the introduction of low-cost depth sensors on the consumer market, several recent methods propose to tackle the problem from a geometric alignment point of view. These generative techniques require to have a geometric template that can be aligned to the sensor data. One of the problems with these methods is that their template is either inaccurate or inefficient. Typically, implicit templates are used because of their ability to efficiently compute distance queries. But they require a lot of primitives to represent a hand accurately, which makes them inefficient. To solve this, [START_REF] Tkach | Sphere-meshes for real-time hand modeling and tracking[END_REF] represent their hand template with a sphere-mesh, a representation that allows good accuracy with very few primitives.

Contributions

With the target of modelling a variety of different objects and shapes with simple primitives, we studied how to combine these primitives. Implicit composition operators are very useful for this but they are complex to design. For this reason, in our first project we introduced a method to design these operators with simple sketches, which is accessible even to non-technical users.

In this first project we often used "tapered capsules", a primitive that has the ability to model a lot of shapes by varying its few degrees of freedom. This primitive is also very efficient for distance queries and, most importantly, it can preserve its volume while being stretched. For these reasons, in our second project we applied this primitive to real-time simulation of rods, muscles and soft-bodies. In our final project, we perform generative modeling in an "indirect" way, by automatically detecting and learning the primitives that occur repeatedly in a given signal. This is done without any supervision regarding the location or shape of the primitives.

In summary, the contributions of this research are:

1. A method to design composition operators. We introduce a method to design functions that combine implicit primitives via simple sketches. This method enables users to apply composition modeling without requiring any technical knowledge. This is achieved by introducing a template that can represent a wide variety of implicit blending operators. This template is optimized to produce the composition described by a user sketch.

2.

A volume-preserving primitive for real-time simulation. We introduce an efficient primitive for volume-preserving deformations. This primitive can be applied to the simulation of rods, muscles and soft-bodies. In order to simulate it in real-time, we extend position-based dynamics by considering scale as a new degree of freedom.

3.

A weakly supervised method to learn and detect primitives. We introduce a deep learning framework to automatically identify the structure of repeating elements in a signal. This is made possible by the introduction of a multi-stage process to train through the non-differentiable top-k operation.

Parts of this research have been published in the following papers: The following chapters detail the contributions of this research. This dissertation is organized so that each chapter contains its own necessary background.

Chapter 2

Composition of Primitives

Abstract

Implicit models can be combined by using composition operators; functions that determine the resulting shape. Recently, gradient-based composition operators have been used to express a variety of behaviours including smooth transitions, sharp edges, contact surfaces, bulging, or any combinations. The problem for designers is that building new operators is a complex task that requires specialized technical knowledge. In this work, we introduce an automatic method for deriving a gradient-based implicit operator from 2D drawings that prototype the intended visual behaviour. To solve this inverse problem, in which a shape defines a function, we introduce a general template for implicit operators. A user's sketch is interpreted as samples in the 3D operator's domain. We fit the template to the samples with a non-rigid registration approach. The process works at interactive rates and can accommodate successive refinements by the user. The final result can be applied to 3D surfaces as well as to 2D shapes. Our method is able to replicate the effect of any blending operator presented in the literature, as well as generating new ones such as non-commutative operators. We demonstrate the usability of our method with examples in font-design, collision-response modeling, implicit skinning, and complex shape design.

Introduction

An implicit representation of a 3D object describes its surface as a set of 3D points on which a scalar function equals a prescribed iso-value [Bloomenthal and Wyvill, 1997]. When modeled or animated, complex objects are defined by assembling their different parts with composition operators, each part being defined by its own scalar function. While the iso-surfaces represent the individual shape of the parts, composition operators control the way they are combined. For instance, the max (min) of two scalar functions produces a union (intersection) operator [START_REF] Sabin | The use of potential surfaces for numerical geometry[END_REF][START_REF] Ricci | A constructive geometry for computer graphics[END_REF] which is the basis of Constructive Solid Geometry (CSG) [START_REF] Requicha | Constructive Solid Geometry. TM[END_REF]; see Figure 2.2. The blending operator, in some cases a simple sum of the combined scalar functions [START_REF] James F Blinn | A generalization of algebraic surface drawing[END_REF], smoothes the sharp transition between parts produced by the union. A core feature of implicit representations is that primitives are combined by simply applying an operator to their respective scalar functions, regardless of their relative positions. This means that no detection or specific treatment for collision is required. This is convenient when the combined primitives are particles of a point-based fluid simulation [START_REF] Ihmsen | SPH Fluids in Computer Graphics[END_REF], or limbs of a character [Vaillant et al., 2013a].

Several composition operators have been proposed, including controlled blending [START_REF] Hoffmann | Automatic surface generation in computer aided design[END_REF][START_REF] Rockwood | The displacement method for implicit blending surfaces in solid models[END_REF][START_REF] Pasko | Function representation in geometric modeling: concepts, implementation and applications[END_REF][START_REF] Hsu | Field functions for blending range controls on soft objects[END_REF], localized blending [START_REF] Pasko | Function-based shape modeling: mathematical framework and specialized language[END_REF], and contact operators that model the contact surface where the combined objects are colliding [START_REF] Cani | An implicit formulation for precise contact modeling between flexible solids[END_REF]; see variety of composition possibilities, they are not commonly used in practice. Some reasons are that meshes are the standard representation for modeling/animation and implicit modeling is not popular on its own, these operators can be computationally intensive, the shape they produce can be unsatisfactory in some cases, and they are often difficult to control by a user.

Recently, gradient-based composition operators [START_REF] Gourmel | A gradient-based implicit blend[END_REF] addressed various unsatisfactory shapes in compositions and computationally expensive operator evaluations. Noting that an implicit surface can approximate a mesh by computing a signed distance field [START_REF] Macedo | Hermite radial basis functions implicits[END_REF], implicit skinning [START_REF] Vaillant | Robust iso-surface tracking for interactive character skinning[END_REF] exploits the automatic contact handling of gradient-based operators on 3D scalar functions for deforming meshes more efficiently when they are animated. This is an example of how implicit modeling/animation tools can be complementary to existing techniques for mesh processing; the current work provides an effective solution to the generalization and intuitive design of free-form gradient-based composition operators.

Composition operators.

The scalar functions f a (x), f b (x) : R 3 → R of two objects can be combined with a binary composition operator g : R 2 → R, and the function f c defines the resulting

object as f c (p) = g(f a (p), f b (p)).
Even though by suitable choices of composition operators g a wide variety of transitions can be obtained, many desired behaviors cannot be captured by any choice of the operator [START_REF] Gourmel | A gradient-based implicit blend[END_REF]. For example, an operator that produces a smooth blend at a transition will also cause a potentially undesired bulging deformation where two objects overlap, as well as premature bulging before the objects make contact; see

f c (p) = g( f a (p), f b (p), θ ) (2.1) where θ = ∠( ∇f a (p), ∇f b (p) ) (2.2)
∠(v, w) being the angle between vectors v and w. In the rare degenerate cases where this angle is undefined, because the gradient vanishes, θ is 0. The domain of the operator g, i.e. the operator-space, is

D = [0, 1] × [0, 1] × [0, π]
This can be understood as using different composition operators for different values of θ. That is, according to whether the two gradients are in opposite, orthogonal, matching directions, or anything in between. Undesired bulging can be resolved, and pre-contact deformations can be disabled, while in the other cases (intermediate angles) the transition can be kept smooth. Four specific problems of implicit modelling (bulging, locality, absorption, topology) were addressed by defining appropriate instances of g [START_REF] Gourmel | A gradient-based implicit blend[END_REF]. The main challenges with gradient-based operators is that designing and fine tuning an operator to obtain some desired effect is a highly technical task, and not possible for non-expert end-users such as 3D artists and modelers. Thus, only predefined and fully parametrized operators can be provided to users and they have to be set in the system by experts.

Contributions.

In this research, we address this usability-gap in composition modeling. We introduce a novel, interactive editing pipeline where the user sketches the desired behavior directly in 2D over one example, and an automatic optimization produces the corresponding operator, transparently to the user. Importantly, the design of an operator and its usage are kept orthogonal:

an operator can be applied in any context (2D or 3D) where the same shape or behavior is desired, irrespective of the example where it was designed. Our method can be applied to any shape for The resulting operator, if applied to the initial pair of primitives, produces the desired behavior and can be applied wherever this effect is desired.

which we can compute a signed distance function. The following contributions were necessary to achieve this result: (1) we introduce a new template to represent operators, which avoids the use of transfer functions such as those proposed in [START_REF] Gourmel | A gradient-based implicit blend[END_REF], and is suitable both in the design and in the application phases;

(2) we present a way to map user-sketches into samples of operator-space D;

(3) we observe that the problem of fitting this template to the samples can be cast as a deformable surface registration problem, and we identify suitable regularizers; (4) we introduce the concept of non-commutative operators, which we show to be useful in certain scenarios, and finally (5) we introduce some novel interesting applications for these operators.

Related Work

Because of the ease of representing arbitrary and changing topologies, CSG, blended and contact surfaces, implicit modeling has some advantages over traditional surface models [START_REF] Marschner | Fundamentals of computer graphics[END_REF]. After presenting the previous works on controllable composition operators and an overview of sketch-based implicit modeling, we review some potential applications of our operators.

On freeform operators. Over the years, several operators have been designed to try to fill the gap between their mathematical formulation and their manipulation by end-users. In aesthetic blends, [START_REF] Alexander | Blending operations for the functionally based constructive geometry[END_REF] optimize the three parameters of an algebraic blending operator to approximate a user's sketch. More complex free-form 2D operators have been defined by blending implicit lines [START_REF] Barthe | Two-dimensional potential fields for advanced implicit modeling operators[END_REF] or by deforming a blending operator [START_REF] Barthe | Controllable binary csg operators for soft objects[END_REF]. These operators are subject to all the limitations solved by gradient-based operators and are not compatible with this operator formulation. Rather than focusing on the operator shape, [START_REF] Pasko | Bounded blending for Function-Based shape modeling[END_REF] and [START_REF] Bernhardt | Implicit blending revisited[END_REF] propose to localize the influence of blending operators on the combined objects by adding an additional 3D scalar function, placed automatically or user-defined. These approaches focus on the definition of where the blending should occur, and no explicit control is performed on the shape of the operator itself. As introduced by [START_REF] Gourmel | A gradient-based implicit blend[END_REF], the shape of the gradient-based operators are defined by a set of 2D profile curves in cylindrical coordinates. This definition is unintuitive to manipulate and restricted in the variety of shapes it can produce. To achieve better skinning with accurate contact deformations, [START_REF] Vaillant | Robust iso-surface tracking for interactive character skinning[END_REF] introduced gradient-based operators discretely computed in a 3D grid, 2D slice by 2D slice, by bi-harmonic interpolation of Dirichlet constraints at boundaries and additional constraints on the iso-value. Their specification relies on a lengthy trial and error process, consisting of editing a set of spline curves for a given θ, as well as a trigonometric transfer function to non-linearly interpolate these curves along the θ axis. Furthermore, an interactive exploration is not feasible, as a re-computation of the computationally intensive fairing optimization is necessary on each update. Finally, their operator construction is tailored to a small set of effects useful in the targeted context (symmetric contacts and bulges), whereas our operators are generic. In our work, we enable the specification of freeform gradient-based operators at interactive rates through the use of 2D annotations, which directly describe the intended user-defined behavior.

Sketch based implicit systems. Sketches have long been recognized as a powerful tool for modeling [START_REF] Igarashi | Teddy: A sketching interface for 3d freeform design[END_REF]. Sketch-based implicit systems added the ability to do blending and CSG with volume models in work such as [START_REF] Singh | Wires: A geometric deformation technique[END_REF][START_REF] De | Blobmaker: Free-form modelling with variational implicit surfaces[END_REF][START_REF] Alexe | Shape modelling by sketching using convolution surfaces[END_REF][START_REF] Tai | Prototype modeling from sketched silhouettes based on convolution surfaces[END_REF], and there have been several examples, including the popular ShapeShop by [START_REF] Schmidt | Shapeshop: Sketch-based solid modeling with the blobtree[END_REF]. [START_REF] Singh | Wires: A geometric deformation technique[END_REF]] shares with us the idea that the final surface shape is modelled by a 3D curve. Closer to our proposal, Karpenko et al.

[2002] built models using an implicit representation based on Radial Basis Functions. Their system used the input stroke to edit a mesh, which would in turn change the implicit representation.

This implicitization approach changes the field locally according to the user's edit. In the above approaches, sketches define the shape in one particular modeling instance. The input sketch in our system defines an operator, which is not tied to the context where it is defined, but can be applied wherever the user desires.

Applications

This work impacts several application domains, offering interesting contributions in each of them.

Character skinning. In the context of character skeleton-driven animations, composition operators have been used to achieve more realistic procedural skin distortions [Vaillant et al., 2013a[START_REF] Vaillant | Robust iso-surface tracking for interactive character skinning[END_REF]. This is also a motivation for our work. With respect to these applications, our approach offers the ability for the designer to intuitively sketch the exact intended deformation (e.g. skin bulging) in one instance, and produce an operator which will reproduce that deformation in real time. This approach is analogous to example-based deformation schemes (see [Jacobson et al., 2014b] or [START_REF] Shi | Examplebased dynamic skinning in real time[END_REF]), but in our case, the exemplar sketch can be drawn in 2D, and the extracted operator can be directly applied to any other joint.

Font design. One potential application of our method is to assist font-design, where glyphs for each letter are the result of composition operators in 2D from a pre-defined skeleton; similarly to [START_REF] Suveeranont | Example-based automatic font generation[END_REF]. The literature on font design is large, and covers many problems which are not addressed here, including automatic construction of the skeletons, or a consistent cross-parameterization for all glyphs, or even a generative manifold of all possible fonts; see [START_REF] Neill | Learning a manifold of fonts[END_REF]. In pipelines where the skeletons of each glyph becomes available, our method offers the ability to control the shape of one glyph (e.g. serifs and joints), and apply it consistently to the entire type-face.

Botanical modeling. In the context of procedural synthesis of botanical models (both realistic and stylized), implicit models have been recognized early as suitable solution, due to their natural ability to recreate smoothly blending branching structures [START_REF] Bloomenthal | Skeletal Design of Natural Forms[END_REF][START_REF] John | Implicit modeling of tree surfaces[END_REF][START_REF] Galbraith | Implicit Visualization and Inverse Modeling of Growing Trees[END_REF]. Using our tool, a user can simply trace the required shape to mimic the geometry of these features, and incorporate this effect into the operator.

Method Overview

A visual outline of our framework is illustrated in Figure 2.5. Our design process begins with the user placing two exemplar implicit primitives in a 2D sketch. The user then annotates the desired blending behavior by sketching a curve. Given this input, an automatic system derives an operator g by solving an inverse optimization problem; see Section 2.6. The resulting operator can then be used both to combine implicit curves in 2D and to combine surfaces in 3D. This observation is crucial, as it allows the user to simply work in 2D to produce operators for 3D modelling. This is even more relevant for the design of contact surfaces, which would otherwise be difficult to edit (or even just to visualize) in 3D, due to self-occlusions.

Feedback loop. In many cases, a single set of user sketches are sufficient to fully determine the operator g that produces the desired result; see Figure 2.5. For more complex cases, an iterative feedback loop can be used to refine the operator and at the same time observe its effect; see Figure 2.6. First, an operator g 0 is constructed from an initial sketch and automatically applied to the exemplar primitives. The resulting 2D drawing provides feedback to the user. The user can then add new sketches, and a new operator is produced from the union of all sketches. This is repeated until a satisfactory shape is returned to the user, limited only by the expressiveness of the gradient-based approach; see Section 2.7. In practice, we found we needed no more than three iterations. For the feedback loop to be interactive, we need to ensure that our algorithms are computationally efficient for the optimization of the operator and for its application.

Background and preliminaries

An implicit model (surface in 3D, or contour, in 2D) is defined by a scalar field-function f , as the set of the points where f assumes a given iso-value. Following the convention from Bloomenthal

and Wyvill [1997], we define the surface as S = {x ∈ R n |f (x) = 0.5} which bounds an interior where f (x) > 0.5. The field-functions of the primitives are, in turn, defined by their skeleton. For example, a sphere is generated by a point skeleton, and a capsule (a cylinder with hemispherical caps) by a line-segment skeleton. Any other shape can be used; see s A central concern in implicit operator design is to ensure smooth blends and avoid normal discontinuities at the boundaries of the support. To this end, where functions such as min or max are used for g [START_REF] Barthe | Two-dimensional potential fields for advanced implicit modeling operators[END_REF], filter fall-off functions are required to be at least C 1 . In our approach, where we fully control the composition function g, this requirement can be completely dropped. Instead, we rely on functions g with built-in smoothness, by defining appropriate value and derivative constraints at the boundaries of its domain D; see Section 2.6.3. This observation allows us to use any monotonic C 0 fall-off function for our primitives. In our examples we opt for a simple linear function controlled by two intuitive parameters: R 1 , the iso-value of the implicit model, and R 2 , the thickness of the support; R 1 is mapped to 1/2 and R 1 + R 2 to 0 (see Figure 2.7, top). In the example of Figure 2.13, the curly branches are obtained by linearly interpolating two values of R 1 along the skeletal curve of the branch.

f a = 0.0 f a = 0.5 f a = 1.0 f b = 0.0 f b = 0.5 f b = 1.0 R 1 R 2 ∇f a ∇f b f b f a θ ab D θ a b f a θ ab s

Intersection and difference.

In this research, we concentrate on union composition operators g, which fuse two primitives into one in some prescribed manner, i.e. g always returns 1 when either of the first two parameters is 1, and 0 when both are 0. Generalization to intersections and differences is straightforward using the same g:

g intersection (a, b, θ) = 1 -g(1 -a, 1 -b, θ) (2.3) g difference (a, b, θ) = 1 -g(1 -a, b, θ) (2.4)
Implicit composition. The creation of complex geometry requires the composition of more than just two input functions. Following the ideas introduced by [START_REF] Wyvill | Extending the csg tree. warping, blending and boolean operations in an implicit surface modeling system[END_REF], we employ different binary operators at each node of a tree, with primitive shapes at its leaves. For example, see Figure 2.1 where different operators are used in cascade.

Capturing user inputs

In our system, a pair of 2D primitives are arranged freely by the user; see Figure 2.6 for an example.

The linear field functions f a and f b of the two primitives are expressed in closed form, and to the user, the two primitives are visualized as closed poly-lines. In Figure 2.7, we also visualize the supports S a and S b of the two shapes. By construction, the target operator can only determine the resulting shape in the area S a ∩ S b , therefore user input strokes are restricted to be inside this area with a stencil mask. The user draws the desired shape of the resulting surface over the input primitives, with one or multiple sketches. For most experiments we employ parametric curves, but any drawing method such as the strokes in Figure 2.20c, can be used. As described below, this is possible as only a sampling of the sketches is required to derive an operator.

Sampling user input.

From the user's sketch, we extract a set of n samples {s 1 , . . . s N }. Each sample s n ∈ R 2 represents a position that the user expects the result/output surface to cross. As illustrated in Figure 2.7, for each sample s n we define a corresponding sample sn in the operator domain D:

sn =      a n b n θ n      =      f a (s n ) f b (s n ) ∠( ∇f a (s n ) , ∇f b (s n ) )      (2.5)
Computations of sn from s n are conveniently fast because functions f a and f b are available in closed form. Their gradient is either available in closed form or approximated by finite differences. The regressed operator g evaluated at sn should return the value (0.5), or in other words g(s n ) = 0.5. Hence, after optimization, the designed blending operator kernel g should interpolate each sample sn .

Sketches over 3D rendering. As a variation, sketches can be drawn over a 3D rendering of the implicit surfaces, seen from arbitrary viewing angle; e.g. see Fig. 2.17. In our prototype, the sketch is in this case assumed to lie on an plane parallel to the image at a predefined depth. In situations where the gradients of the primitives are not co-planar, there is no ideal plane to sketch on and results could be unintuitive. In this case, it is often better for the user to sketch on a simpler configuration of primitives and to then apply the operator to the original scene. Our framework requires no further adaptation to deal with this case. A more advanced interface could identify depth automatically, for instance by projecting the first point of the sketch on the shape, similarly to the sketch-based modeling interface proposed by [START_REF] Bernhardt | Matisse: Painting 2d regions for modeling free-form shapes[END_REF].

Fitting the composition operator

The blending operator requires a 3D function g to be defined over its entire space D, however samples collected from user's strokes only define the behavior of g in a small subset of D. On the other hand, many characteristics of the function g are known a-priori, such as its general shape (Sec. 2.6.1), its continuity requirements (Sec. 2.6.2) and the values on the boundaries of D (Sec. 2.6.3). We approach this reconstruction problem in two steps: first, we identify the set of 0.5 values of g, as a parametric surface P embedded in D, which is fitted to the samples (Sec. 2.6.4); we then compute a 3D lattice covering the domain of g by propagating the iso-values in P (Sec. 2.6.5).

The final operator is then evaluated by tri-linear interpolation of the lattice values.

The operator template

In previous work [START_REF] Gourmel | A gradient-based implicit blend[END_REF]Vaillant et al., 2013a[START_REF] Vaillant | Robust iso-surface tracking for interactive character skinning[END_REF], g(a, b, θ) is formulated as a collection of two dimensional functions for a few particular values of θ, each independently defined as a curve defining the portions of the domain where g evaluates to (0.5). These cases can be interpreted as axis aligned slices of the domain D. In our work, we conveniently represent the operator g as one surface P embedded in D representing its (0.5) iso-value; see Figure 2.5b. This approach allows us to define a template for the surface P , designed to represent a wide class of useful operators: sharp unions (see Figure 2.8a), smooth blends (see Figure 2.8b), articulated contact (see Parametric representation. Our template for P is a surface made by a pair of third-order B-Spline patches, each with I × J control points, joined at their boundary and arranged as in Figure 2.8a. The surface P is fully determined by the positions of the control points p a i,j and p b i,j in D. We found I=J=5, for a total 45 distinct control points, to provide the necessary expressiveness while avoiding excessive redundancy. The continuity of P at the junction is enforced by imposing ∀i ∈ [1..I] : p a i,J = p b i,J , while other boundary and regularization constraints will be discussed later.

Template expressiveness

A fundamental characteristic of our operator template lies in its expressiveness. In particular, in Figure 2.21 we illustrate how, to the best of our knowledge, all results obtained by any composition operators which have been proposed in the literature can be expressed by our template. We now detail how several operators can be realized by properly deforming our template.

Sharp creases: if a sharp crease is desired in the transition between the two operand surfaces (e.g. with union) then g needs to break C 1 continuity, and consequently P must have a normal discontinuity. In our template, this discontinuity is easily accommodated by construction at the junctions between the two splines.

Smooth blends: if g is required to generate surfaces without any sharp creases in the transition between the two operand surfaces, it needs to be C 1 continuous, and consequently, surface P must be smooth, including at the junction between splines. This can be easily obtained by

aligning the points {p a i,J-1 , p a i,J = p b i,J , p b i,J-1 }.
Contact surfaces: another important feature of operators is the ability to produce contact surfaces at the transitions between the two primitives [START_REF] Cani | An implicit formulation for precise contact modeling between flexible solids[END_REF][START_REF] Vaillant | Robust iso-surface tracking for interactive character skinning[END_REF]; see ). Clearly, non-commutative operators, such as those in Figure 2.8d can also be represented.

In summary, our template can recreate each of the features above. If required, our system allows users to explicitly enforce these conditions; however, it is not strictly necessary to do so, as surface P will naturally conform to these conditions whenever they are suggested by the user data as a result of the fitting process. This observation can aid the design of user interfaces based on our method.

Boundary conditions

When f a (p) = 0, the point p is outside the (compact) support of f a , and hence beyond its range

∀θ, ∀a : g(a, 0, θ) = a and g(a, 1, θ) = 1 ∀θ, ∀b : g(0, b, θ) = b and g(1, b, θ) = 1 (2.6)
Further, to achieve (normal) shading smoothness at the boundaries of the supports, we must enforce C 1 continuity of the blending operation by vanishing the derivatives:

∀θ, ∀a : ∂g ∂b (a, 0, θ) = 0, ∀θ, ∀b : ∂g ∂a (0, b, θ) = 0 (2.7)
As θ represents the unsigned angle between the two gradient vectors, the function g is implicitly mirrored at both ends [0, π] of its third parameter. Therefore, to achieve C 1 continuity we also impose:

∀a, ∀b : ∂g ∂a (a, b, 0) = 0 ∀a, ∀b : ∂g ∂b (a, b, π) = 0 (2.8)

The constraints above translate into Dirichlet and Neumann boundary conditions for surface P , which can be conveniently expressed in terms of linear hard constraints on its control points.

To enforce Equation 2.6, the control points on opposite ends of P are constrained to lie on the two line segments at the boundary of D: (0.5, 0, θ) and (0, 0.5, θ). To enforce Equation 2.7 and Equation 2.8, the control points neighboring the boundaries of P must be constrained to be vertically/horizontally aligned to the boundary control points.

Surface registration

In this phase, we fit the template parametric surface P to the collected samples si . P maps each point uv = (u, v) of its 2D parametric space into a position P (uv) ∈ D, as determined by the control points {p k i,j }. We start with an initial guess, which corresponds to a union operator: a surface P where the two B-spline patches are simply planar and reciprocally orthogonal; see Figure 2.8a. We then non-rigidly register the template onto the samples following the approach in [START_REF] Bouaziz | Modern techniques and applications for real-time non-rigid registration[END_REF]; see Figure 2.9. This optimization consists of the alternation of two

t = 1 t = 3 t = 4 t = 6 t = 9
Figure 2.9: We illustrate the sketches in operator space, and a few iterations of the operator registration optimization in Equation 2.10. The manifold surface folds onto itself to be able to reproduce a non-manifold configuration. At the same time, the way we enforce the symmetry constraint helps the optimization to avoid undesirable self-intersections. In this example, the optimization converged after 9 (two-steps) iterations; none of our experiemnts required more than 20.

optimization steps:

local: arg min (uvn) P (uv n ) -sn 2 , ∀s n (2.9)
global: arg min

{p k i,j } E match ({p k i,j }, {uv n }) + E priors ({p k i,j }) (2.10)
That is, our local-global ICP optimization first computes closest-point projections, and then it globally modifies the surface control points, resulting in a non-rigid deformation of our surface. For sample points corresponding to a contact surface, both sides of the template must be separately fitted to these. These sample points could be labeled automatically but we let the user provide this information. To efficiently implement the local step, we first triangulate P (we employ a resolution of 40 × 40), and use a regular volumetric grid to accelerate closest-point queries from sn to the triangles. We also constrain each of the I rows of control points to lie at a constant equally spaced θ values. The global step is implemented as one Least Squares minimization, as both energy terms are quadratic in the variables; E match represents the data-to-model error, while E prior accounts for shape-priors as well as optimization regularization. Further implementation details are available in our publicly released source code.

Matching term. The data-to-model error is computed as the averaged squared distances of samples from the tangent planes of their projections onto P :

E match = 1 N n [n n • (s n -P (uv n ))] 2 (2.11)
where n n is the normal at P (uv n ). This point-to-plane metric leads to better convergence compared to point-to-point errors. The term is a quadratic function of the variables, because P (uv n ) and n n are both constants in the global step.

Priors. The prior energy includes several terms weighted by parameters. These terms enforce:

operator fairness, potential contact constraints, and regularization of the optimization:

E prior = w fair E fair + w contact E contact + w tikh E tikh (2.12)
Fairness prior. We penalize oscillations in the surface with a bi-harmonic energy defined on the control points p k i,j :

E fair = i,j ∆ uv p k i,j 2 (2.13)
where ∆ uv represents the 2D Laplacian operator defined in the uv domain, adapted to have control vertices stored in matrix form p k i,j . This regularizer has multiple advantages: (1) in under-sampled areas it ensures a smooth interpolation, (2) it prevents over-fitting and, (3) as the sketches only provide a sparse sampling in D, it regularizes our optimization ensuring the problem remains well-conditioned. Additionally, these fairness energies are known to penalize surface fold-overs [START_REF] Botsch | A remeshing approach to multiresolution modeling[END_REF]. Following [START_REF] Li | Global correspondence optimization for non-rigid registration of depth scans[END_REF], we start by a strong enforcement of fairness to avoid local minima, and progressively relax this constraint to allow the surface to eventually closely fit to the data. Specifically, we employ the weight scheduling

w f air = 10 3 • 2 -t + 10 -4
, where t is the iteration number.

Full-contact prior. When we detect that the user sketch corresponds to a full-contact (e.g. see Figure 2.21e, where two objects are fully separated by a contact interface), we also enable a prior that ensures the seam control points connecting the two patches of our template project on the boundary of the domain D:

E contact = (i,j)∈S {nm} n i • (p i,j -[1, 1, 1]) 2 n 1 = [1, 0, 0] n 2 = [0, 1, 0] (2.14)
Because of the multiplication , this energy is non-linear, hence when computing its gradient we only enable the term that has the smallest point-to-plane residual.

Tikhonov regularization.

As the fitting energy is linearized within each local step, we avoid overshooting with a mild Tickhonov regularizer, by setting w tikh = 10 -3 which penalizes displacements from the previous solution:

E tikh = i,j,k p k i,j (t) -p k i,j (t -1) 2 (2.15)

3D-Lattice filling

Once the surface P describing the (0.5) iso-values of g is defined, the next step is to regularly sample its domain D and assign scalar values on each cell of the lattice. This is executed in three sub-steps:

Step 1 -Initialization: we assign the voxels immediately surrounding P with the signed distance from either of its two patches, by rasterizing them over the 3D lattice. Grid receiving two values are set to the greatest (i.e. most internal) one; this ensures robust evaluation when the two patches are coplanar, as in the case of contact operators, as the contact surface is in the interior of the object.

Step 2 -Boundary assignment: for each θ slice of D, we assign the values to the four sides et al., 2013], as well as the (bottom) implicit-contact from [Vaillant et al., 2013a]. {(0, * , θ), (1, * , θ), ( * , 0, θ), ( * , 1, θ)} according to the boundary constraints in Equation 2.6.

k = 0 . 2 k = 1 . 2 k = 0 . 4 k = 0 . 6 k = 0 . 8 k = 1 . 4 k = 1 . 6 k = 1 . 8 k = 0 . 0 k = 1 . 0 g a g b
Step 3 -Value propagation: the values assigned in the previous two sub-steps are diffused over the remaining portions of D by solving a bi-harmonic fairing optimization (where the values set in previous steps act as hard constraints).

Step 3* -Efficient value propagation:

Step 3 is time-consuming, as its complexity is cubic in the linear resolution of the lattice. Fortunately, this is only necessary when the operator must be used in successive compositions. When only two objects need to be combined, as during the operator design process, only the lattice values surrounding the (0.5) surface are relevant, and the other ones can be safely disregarded. This observation drastically reduces the latency of the feedback loop, effectively enabling the interactive design discussed in Section 2.3.

Evaluation

We evaluate our work by verifying the expressiveness of our sketched operators, the controllability of results, as well as the possibility of interpolating between different operators.

Expressiveness. We empirically validate our approach by demonstrating a variety of effects. We collect all the fundamental types of implicit blending operator that have been proposed in the literature (to the best of our knowledge), and verify that our template can be optimized to express the same behavior. To this extent, Figure 2.21 reports a number of representative images from the literature, and the 2D sketch necessary to generate the desired operator. We visualize the optimized (deformed template) operator, and the result of its application on the 2D input geometry.

We also show how all our custom operators extend to 3D in a completely straightforward fashion.

Our sketches can be used to: (a) represent traditional CSG operations [START_REF] Sabin | The use of potential surfaces for numerical geometry[END_REF]. (b) smoothly blend two primitives [START_REF] James F Blinn | A generalization of algebraic surface drawing[END_REF] and [START_REF] Ricci | A constructive geometry for computer graphics[END_REF], (c) perform bulge-free blend [START_REF] Gourmel | A gradient-based implicit blend[END_REF]. (d) avoid premature-blending [START_REF] Gourmel | A gradient-based implicit blend[END_REF], (e) model bulge-on-contact between two objects [START_REF] Cani | An implicit formulation for precise contact modeling between flexible solids[END_REF], as well as (f,g) represent the partial-contact from implicit-skinning [Vaillant et al., 2013a]. In fact, the expressiveness of our operators go beyond the capabilities of those proposed in the literature, and allows us to extend the contact operators pioneered by [START_REF] Cani | An implicit formulation for precise contact modeling between flexible solids[END_REF],

towards the representation of asymmetric-contact, without having to edit the input scalar fields; see Figure 2.10.

Controllability and interpolation.

As our algorithm receives user sketches as input, the behavior of the operator is easily controllable. In Figure 2.11, we demonstrate how several variants of blending and contact operators can be faithfully reproduced by specifying the desired behavior with a simple 2D sketch. Another form of composition control can be obtained by blending two distinct composition operators. As illustrated in Figure 2.8, the topology of the template is consistent across all of our examples. This allows us to blend operators through simple linear interpolation/extrapolation of its control points; see Figure 2.12. no effect mercury water cross-section Figure 2.15: Applying implicit contact operators to the output of a fluid simulation produces a noticeable improvement in visual quality. This effect is efficient to apply, as many fluid solvers represent liquids with implicit functions. While our operator can apply to any geometry, we choose a simple and familiar scene that can be more easily appreciated.

R (1) 1 R (2) 1

Applications

We present several example applications that benefit from the sketch-based generation of blending operators. The input sketches can either be drawn freely, or rotoscoped over annotated images.

Botanical modeling. Our curly tree shown in Figure 2.13-top is inspired by the images returned by the search query "curly tree". The placement and shape of primitives can be produced by any procedural grammar-based algorithm [START_REF] Lienhard | In Design Transformations for Rule-based Procedural Modeling[END_REF]. Our technique allows for the automatic generation of controllable smooth-blends connecting branches to each other. Two user strokes are sufficient to define the desired operator behavior, and the result is applied to a large quantity of compositions. In this case a single operator was used, but variability can be easily achieved by interpolating multiple variants (see Figure 2.12) with random weights; a 3D example of this process is visualized in Figure 2.17. In Figure 2.13-bottom, we apply our Figure 2.17: (a) A procedural tree generated by [START_REF] Lienhard | In Design Transformations for Rule-based Procedural Modeling[END_REF] as a collection of spheremeshes [START_REF] Tkach | Sphere-meshes for real-time hand modeling and tracking[END_REF]. (b) a traditional blending operator results in unwanted bulges (e.g. Blender's metaballs), (c) our blending operator designed via a simple 2D sketch is propagated to all branches in the hierarchy.

sketch-driven operators to the procedural generation of botanical leaves. The skeletal structure of the leaf is generated via an appropriate grammar, as in the tree example. In this setup we experiment with two different sketches. In Figure 2.13d, the sketch is the entire boundary of the leaf from Figure 2.13c, while in Figure 2.13f the user has sketched the operator himself. The optimized operators result in blended geometry closely mimicking the geometry of natural leaves (we superimpose the original texture to make this more apparent).

Approximating contact behavior. In Figure 2.14, we approximate the behavior of a rubber ball in motion through the design of two operators. The first distorts the shape of the ball, modeled via a "ghost" primitive whose position is determined according to the velocity vector.

The second operator captures the contact behavior between the projectile and the deformable wall. In our stop-motion illustration the distortion was interpolated, while the contact behavior was extrapolated. Asymmetric contact operators can also be used to approximate the fine-scale behavior of liquids when they come in contact with surfaces coated with different materials; see Figure 2.16. Obtaining such effects through fluid dynamic simulation is difficult as it requires a volume-preserving solver and careful implementation of boundary conditions; [START_REF] Wang | Water drops on surfaces[END_REF]. Such effects, due primarily to surface tension, can also be approximated by evaluating a curvature flow on surfaces [START_REF] Thürey | A multiscale approach to mesh-based surface tension flows[END_REF]. Our approach approximates these effects without the need for complex physical simulation, and allows their application in a lightweight fashion as a post-processing step; see Figure 2.15. Implicit vector font design. Another application is the creation of font libraries, in particular as our method allows a user to explore the design space interactively; see Figure 2.18. In this example, we use multiple sketches on the letter T and A to derive one operator, reproducing the desired behavior for orthogonal and non-orthogonal configurations. Again, the composition operators is derived from a few exemplars and then applied to all other cases. The user can explore different variants for a character, and export the blending operator to create a self-consistent font library. An example showcasing the extensibility of font design to three dimensions is illustrated in Figure 2.1. In typography, another common task is the application of serifs to characters. In this example, we apply our composition operators to this task. To achieve this, the yellow primitive removes ink from the pink primitive, with a difference operation (see Equation 2.4). This setup provides the user with a fast and easy method to prototype different serif styles, as illustrated.

Because our operator kernel P is an algebraic surface, the operator g can be generated at any scale; hence, as long as the input primitives can also be expressed in algebraic form (as in the examples shown), our fonts can be synthesized at any desired resolution.

Conclusions

Implicit models have recently re-emerged as a powerful modeling technique with the introduction of gradient-based operators [START_REF] Gourmel | A gradient-based implicit blend[END_REF], solving several problems that were often regarded as intrinsic to implicit modelling. While more expressive the new operators are difficult to control.

In this work we have introduced a solution based on an inverse approach which allows sketch-based design of new operators. The designer directly describes the desired effect of the operator on the resulting surface, without having to understand the mathematics of the operator space. We have shown a number of practical applications demonstrating that our approach enables a designer to quickly take advantage of the powerful nature of implicit modelling. In addition, our approach also serves as a conceptual tool to investigate the space of possible gradient-based operators. Using it, we were able to create a new class of non-commutative operators; see Figure 2.10.

Limitations and future work. The expressiveness of the gradient-based implicit operators, defined as in Equation 2.1, while being superior to the traditional 2D operators, is still insufficient for some cases. For example, regardless of how it is designed, no such operator can reproduce the shape sketched in fitting process provides an "average" shape that can be unsatisfactory when the input sketches are too different. In the future we intend to further generalize blending operators, to augment the range of achievable effects; e.g. anisotropic operators. We have shown that operators based on a 3D domain can be designed with 2D sketches; the same approach could be applied to higher dimensional operators, which have additional parameters.

29 The applicability of our custom operators is not limitless either. While they can be applied to any implicit model, the fall-off functions are assumed to be roughly similar, otherwise the shape produced by the operator can diverge arbitrarily from the user's sketch; e.g. if one of the two support radii R 2 is vastly different.

C 1 C 0 (a) (b) (c)
Another inherited limitation of gradient-based operators is that they can generate artifacts where scalar fields exhibit gradient discontinuities. This rarely occurs in practice.

A different issue is that our operators, in general, lack associativity. Associative binary operators are desirable because they constitute a natural definition of operators working on any number of operands. An interesting open question is whether the input primitives could be inferred along with the operators. Finally, it would be interesting to handle a sketch on more than two overlapping primitives.

(a) union [START_REF] Sabin | The use of potential surfaces for numerical geometry[END_REF] (b) smooth blend [START_REF] James F Blinn | A generalization of algebraic surface drawing[END_REF] (c) bulge-free blend [START_REF] Gourmel | A gradient-based implicit blend[END_REF] (d) blend+union [START_REF] Gourmel | A gradient-based implicit blend[END_REF] (e) contact+bulge [START_REF] Cani | An implicit formulation for precise contact modeling between flexible solids[END_REF] (f) skinning [START_REF] Vaillant | Robust iso-surface tracking for interactive character skinning[END_REF] (g) skinning+bulge [START_REF] Vaillant | Robust iso-surface tracking for interactive character skinning[END_REF] target image primitives+sketches optimized operator applied in 2D applied in 3D

Figure 2.21: Our sketch-based construction method can be used to effortlessly produce every gradient blend operator proposed in past literature (pioneering citation is indicated). Starting from images depicting the final results of each operator (leftmost column), we simply sketched over 2D contour (second column). From these sketches, our system reverse-engineers the operator (central column). This can be applied (in 2D or 3D alike, rightmost columns) to reproduce the same results. Our target images are courtesy of [START_REF] Gourmel | A gradient-based implicit blend[END_REF] and [START_REF] Vaillant | Robust iso-surface tracking for interactive character skinning[END_REF].

31

volumetric deformation. The source code of our physics engine will be openly available upon publication.

Introduction

In Rods with volume. Within this technological landscape, of particular relevance to our work is the modeling of elastic "rods" [START_REF] Spillmann | Corde: Cosserat rod elements for the dynamic simulation of one-dimensional elastic objects[END_REF][START_REF] Umetani | Position-based elastic rods[END_REF][START_REF] Kugelstadt | Position and orientation based cosserat rods[END_REF]. These models extend "ropes" by augmenting each segment composing the rod with an orthogonal coordinate frame, hence allowing the modeling of torsion on top of stretching/bending. Our VIPER rods extend these formulations by accounting for volume preservation, a phenomena not modeled by existing position-based rod models. Many interesting phenomena require this constraint (e.g. soft-bodies, fluids). For example, water is the largest constituent of most animal tissues (≈ 80% in muscles) hence modeling quasi-incompressible phenomena is of critical importance to achieve believable motion. We address this problem by adding a per-vertex scaling degree of freedom -a measure of the local rod cross-sectionand optimizing for this quantity within the physics solve. Our rod segments are hybrid surface representations, they are explicitly parameterized by the position and scale/radius of their vertices, but their surface is defined implicitly. This hybrid structure makes them particularly well suited for efficient collision detection/resolution [START_REF] Green | Particle simulation using CUDA[END_REF].

Anatomical modeling. Such a physical model not only satisfies our fixation in efficiently simulating rubber bands, but has immediate applications towards the modeling of muscles. As illustrated in Figure 3.2, striated skeletal muscles1 in human bodies can be represented as a collection of fibers surrounded by connective tissue (i.e. fasciae). Simulating these muscle types efficiently is a fundamental problem, as they represent from 36% to 42% of the average human body mass. In this paper, we propose to efficiently model muscles as a structured collection of volume-preserving rods. This new model can also be interpreted as a generalization of the static volumetric primitives in Implicit Skinning [Vaillant et al., 2013b], where skin can then be efficiently modeled as a triangular mesh sliding on an implicit iso-surface defined by our fibers. The VIPER primitive is designed to be integrated with other existing components to produce a complete character representation. These include representations of the skeleton, fat, skin, etc. Volumetric simulation. In the industry, volume-preserving simulation is typically performed by discretizing the interior of the object with tetrahedra or using an approximating cage. However, to the best of our knowledge, even modestly sized models require a lengthy preprocessing (e.g. a large scale eigen-decomposition for computing the deformation modes; see [START_REF] Barbič | Real-time subspace integration for st. venant-kirchhoff deformable models[END_REF]) before real-time simulation becomes possible. For example, while visually striking, computing the simulation in Figure 3.3 requires a two-pass optimization ( 1 muscle, 2 skin). 1 Muscles are discretized with 51k tets sharing 21k vertices, and each of their 4 steps of simulation requires 3.2 seconds. 2 Skin (78k vertices) is simulated as cloth layered on top of fat (having 68k tets sharing 18k vertices), where each of the necessary 4 substeps of simulation requires ≈ 35 seconds of compute. Overall, this cumulates to ≈ 2.5 minutes of compute/frame. Offline simulation can be exploited to learn sub-spaces, which then enables dynamic deformations in real-time; see [START_REF] Xu | Pose-space subspace dynamics[END_REF]. However, once learnt, the dynamic behavior of the model is "baked". Clearly, this is an obstacle towards the ultimate goal of truly interactive physics simulation, and, consequently, interactive modeling. The VIPER primitive not only allows the real-time volumetric simulation of complex anatomical structures, but also provides a viable alternative to cages as a compact control structure for soft-body deformations.

Automatic fiber-bundle modeling. Rod-based representations of muscle fascia are not commonly available -typical asset databases contain tetrahedral mesh models instead. While artists sometimes model main characters to the level of interior muscles, this effort is expensive and not justified for background characters. Thus we also introduce a technique to convert existing assets with minimal user intervention. Our solution builds fiber bundles by first creating a set of slices through the muscle then performing an iterative optimization to interpolate these slices with a given number of rods such that they approximate the input surface well.

Contributions.

Our fundamental contribution is the design of a novel real-time physics engine for soft-body dynamics. Our system presents several sub-contributions:

• We enrich position-based solvers by introducing a new volume-preserving cosserat rods model and associated constraints.

• We demonstrate how these primitives, when assembled into fiber-bundles, are effective in efficiently modeling muscles.

• We introduce a technique for conveniently creating fiber-bundles models from existing simulation assets.

• We introduce the use of VIPER rods as efficient deformation proxies for soft-body deformation.

Related Work

We overview the literature from different angles. We recap example-based modeling frameworks that are commonly used in digital production, as well as recent efforts towards the use of simulated anthropomorphic models. We also review methods that attempt to "emulate" them via geometric processes, and finally processes to calibrate a given model to a target. with no consideration of underlying volumetric structures, often resulting in non-physically realistic effects such as the candy wrapper problem of (LBS) Linear Blend Skinning [Jacobson et al., 2014a].

While these artifacts can be resolved [START_REF] Kavan | Skinning with dual quaternions[END_REF][START_REF] Binh | Real-time skeletal skinning with optimized centers of rotation[END_REF], skinning solutions lack details such as tendons, muscle bulges, wrinkles, and volume preservation. Example-based approaches such as (PSD) Pose-Space Deformation [START_REF] Lewis | Pose space deformation: A unified approach to shape interpolation and skeleton-driven deformation[END_REF][START_REF] Kurihara | Modeling deformable human hands from medical images[END_REF] and BlendShapes [START_REF] Lewis | Star: Practice and theory of blendshape facial models[END_REF] interpolate artist-sculpted shapes to emulate all these effects. However, "producing effects such as skin sliding over underlying structures, or the collision effects visible parts of the body press together, can require considerable skill and weeks to months of sculpting depending on the required quality" [START_REF] Yuen | Personal communication. Head of Creatures, Method Studios[END_REF]. Such data-driven models can be learnt from measurements for both static [START_REF] Loper | Smpl: A skinned multi-person linear model[END_REF], as well as dynamic [START_REF] Pons-Moll | Dyna: A model of dynamic human shape in motion[END_REF] humans, but they hardly generalize outside of their corresponding training domains.

Physically-based anthropomorphic models. Physically-based simulation of characters has a long history [START_REF] Terzopoulos | Physically-based facial modeling, analysis, and animation[END_REF][START_REF] Scheepers | Anatomy-based modeling of the human musculature[END_REF][START_REF] Sifakis | Automatic determination of facial muscle activations from sparse motion capture marker data[END_REF] but, due to its high computational cost, it has only recently began to see practical use. For skeletal muscle deformation, [START_REF] Lee | A survey of modeling and simulation of skeletal muscle[END_REF] provides an excellent overview of the field, in regards to which [START_REF] Saito | Computational bodybuilding: Anatomicallybased modeling of human bodies[END_REF], with its ability to reach near-interactive runtime, can be considered the state-of-the-art. Similarly to blendshape generation, training data can be exploited to generate efficient low-dimensional simulations [START_REF] Xu | Pose-space subspace dynamics[END_REF][START_REF] Schumacher | Efficient simulation of example-based materials[END_REF][START_REF] Botsch | Projective dynamics: Fusing constraint projections for fast simulation[END_REF], enabling physically-based digital characters in production settings [START_REF] Clutterbuck | A physically based approach to virtual character deformations[END_REF]Ziva Dynamics]. A shortcoming of these methods is the requirement that the training set encompasses samples of all configurations of the object/character that that will be needed.

Physically based approaches also permit a decomposition into layers -skin, muscle, fat, and bones -enabling appropriate algorithms to be used for each. Highly relevant to our work is the simulation of skin layered over volumetric primitives pioneered by [START_REF] Li | Thin skin elastodynamics[END_REF], and its realization in commercial software [Vital Mechanics], as well as other existing research [START_REF] Saito | Efficient and robust skin slide simulation[END_REF], industry [Ziva Dynamics], and proprietary solutions [START_REF] Clutterbuck | A physically based approach to virtual character deformations[END_REF] to this complementary problem. Recently, [START_REF] Romeo | Muscle Simulation with Extended Position Based Dynamics[END_REF] proposed the use of PBD for muscle simulation, but its ≈ 40s/frame of processing time makes it unsuitable to interactive applications.

Non physically-based anthropomorphic models. Recent efforts have been made to create alternative representations for sub-skin volumetric models (i.e. representations of muscle, fat, and bones). For example, Maya Muscle [START_REF] Comet | Maya muscle[END_REF] represents muscles via NURBS that drive the skin via LBS, whose volume is artist-driven in a PSD fashion. However, due to the explicit nature of NURBS, collisions are expensive, hence the performance of the system does not scale well in the complexity of the model. Rather than driving skin via LBS, Implicit Skinning [Vaillant et al., 2013b] lets it slide on top of implicit surfaces via optimization. These surfaces are defined by blending components that abstract entire body parts (i.e. union of bone, muscle, and fat). In contrast to our work, note that this model is purely kinematic (i.e. no physics). Another relevant class of methods simplifies anthropomorphic components even further. For example representing an entire arm as two tapered capsules is advantageous for arm vs. cloth collision detection [START_REF] Muller | NVIDIA PhysX SDK 3.4.0 documentation[END_REF]. Sphere-Meshes generalize these representations, and have recently been used to approximate geometry [START_REF] Thiery | Sphere-meshes: shape approximation using spherical quadric error metrics[END_REF], and track its movement in real-time [START_REF] Tkach | Sphere-meshes for real-time hand modeling and tracking[END_REF].

Calibrating anthropomorphic (volumetric) models. [START_REF] Hamadi | Anatomy transfer[END_REF] pioneered the transfer of anatomical structures from a template to a target human via approximate deformation models of soft tissues, and [START_REF] Zhu | Adaptable anatomical models for realistic bone motion reconstruction[END_REF] calibrated these models from a set of RGBD images capturing a human in motion. Analogously to these method, physically inspired models [START_REF] Saito | Computational bodybuilding: Anatomicallybased modeling of human bodies[END_REF] can be calibrated to a set of 3D surface scans [START_REF] Kadleček | Reconstructing personalized anatomical models for physics-based body animation[END_REF]. Of particular relevance to our method is the fiber estimation technique pioneered by [START_REF] Hon | Skeletal muscle fascicle arrangements can be reconstructed using a laplacian vector field simulation[END_REF] employed in [START_REF] Saito | Computational bodybuilding: Anatomicallybased modeling of human bodies[END_REF]. While [START_REF] Saito | Computational bodybuilding: Anatomicallybased modeling of human bodies[END_REF] is interested in deriving the anisotropic deformation frame, we require an explicit decomposition of the muscle in fiber bundles.

Generalized Rod Parameterization

Our physical model of Cosserat rods consists of a smooth parametric curve in 3D space c(z) : p(x, y, z) ≡ p(q, z) = c(z) + s(z) R(z)q (3.1)

[z 0 , z 1 ] → R 3 . An orthogonal frame D(z) ∈ R 3×3 is
where q = [x, y, 0] T is a point on a disc D(z) of radius r(z), aligned with the xy plane, and centered at its center of mass. Similarly, a second function maps the rod from parameterization to its deformed configuration:

p(x, y, z) ≡ p(q, z) = c(z) + s(z)R(z)q (3.2)

Variational Implicit Euler Solver

Our solver is based on the variational form of implicit Euler integration [START_REF] Martin | Example-based elastic materials[END_REF]. The physical model evolves through a discrete set of time samples, with simulation step size h. At time t the deformed volume is defined as p t (x, y, z) and the velocity as ṗt (x, y, z). The rest pose is defined as p(x, y, z) = p 0 (x, y, z). The mass m is assumed to be uniform over the rod. The sum of the external forces is denoted as f ext (x, y, z). We will now drop the indexing (x, y, z) whenever possible to improve readability. We consider position dependent internal forces such that the sum of the internal forces is (3.3) where W i (p, p) is a potential energy function, and K i is a stiffness matrix which we assume to be uniform over the rod, and the notation x 2 A means x T Ax. We can then write implicit Euler as an optimization describing the compromise between an inertia potential and the elastic potentials:

f int (x, y, z) = -1 2 i ∇ p W i (p, p) 2 Ki ,
min {ct,st,Rt} z1 z0 D(z) m 2h 2 p t -pt 2 2 inertia + 1 2 i W i (p t , p) 2
Ki elastic dx dy dz (3.4) With pt we indicate the inertial prediction for p t , i.e., its next position in absence of internal forces:

pt = p t-1 + h ṗt-1 + h 2 m f ext , (3.5)
where ṗ(q, z) = ċ(z) + ṡ(z)R(z) + s(z) Ṙ(z) q.

Discretization. We discretize the curve in the parametrized domain using a set of m + 1 points

{z [0] , . . . , z [m] } connected using m piecewise linear elements of length {l [1] , . . . , l [m] }; see Figure 3.4.
We can approximate the curve integral by integrating over these piecewise linear elements using the midpoint rule. For the integration we also define a set of m midpoints {z [.5] , . . . , z [m-.5] }. A point on a midpoint cross section is then parametrized as:

p(q, z [j-.5] ) = c(z [j-.5] ) + s(z [j-.5] )R(z [j-.5] )q (3.6) c(z [j-.5] ) ≡ 1 2 c(z [j-1] ) + c(z [j] ) (3.7) s(z [j-.5] ) ≡ 1 2 s(z [j-1] ) + s(z [j] ) (3.8)
This is similar to the staggered grid discretization of previous work [START_REF] Spillmann | Corde: Cosserat rod elements for the dynamic simulation of one-dimensional elastic objects[END_REF][START_REF] Grégoire | Interactive simulation of one-dimensional flexible parts[END_REF], where the frames R are stored at the midpoints. Contrary to previous approaches, our model also has a scale that we store at the endpoints of the linear elements. We can now rewrite Equation 3.4 as:

min {ct,st,Rt} m j=1 l [j] D(z [j-.5] ) m 2h 2 p t -pt 2 2 + 1 2 i W i (p t , p) 2
Ki dx dy (3.9)

Elastic Potentials

In this section we detail the elastic potentials used to simulate our volume preserving rods.

Strain potential. We define the strain at a midpoint as

E strain (z [j-.5] ) = D(z [j-.5] )
R T ∇p -RT ∇p 2 Ks dx dy, (3.10) where K s = [ k x s e x k y s e y k z s e z ] is a diagonal stiffness matrix and [e x e y e z ] is the standard basis. ∇p and ∇p denote the deformation gradients, i.e., the Jacobian matrices of the deformation functions [START_REF] Sifakis | Fem simulation of 3d deformable solids: A practitioner's guide to theory, discretization and model reduction[END_REF]. As p and p map R 3 to R 3 , the Jacobian matrices are 3 × 3.

The Jacobians are not rotational invariant so we rotate them back to the parametrization domain to be able to compare them on a common ground. As explained in Appendix A.1 integrating the strain energy (3.10) leads to (3.14) where H s = [ πr 4 k z s e x πr 4 k z s e y πr 4 (k x s +k y s )e z ] is the second moment of area of a disc scaled by the stiffness, and the Darboux vector [START_REF] Kugelstadt | Position and orientation based cosserat rods[END_REF]

E strain (z [j-.5] ) = πr 2 k z s ∇ z c -w wT ∇ z c 2 2 (3.11) + πr 2 (k x s + k y s )(s -s) 2 (3.12) + 1 4 πr 4 (k x s + k y s )(∇ z s -∇ z s) 2 (3.13) + sΩ -s Ω 2 Hs ,
is denoted by Ω = [Ω u , Ω v , Ω w ] T .
Note that we retrieved similar energies in previous works augmented by our additional scale degree of freedom. The energies (3.11) and (3.12) respectively measure the stretch along the curve and the cross section, while (3.13) measures the variation of scale across sections, and (3.14) measures bending/twisting. Interestingly, Equation 3.13 can also be interpreted as a measure of surface stretch.

We use an additional energy measuring the second order variation of scale complementing (3.13) with a measure of surface bending

E bending (z [j-.5] ) = 1 4 πr 4 (k x b + k y b ) ∇ 2 z s -∇ 2 z s 2 (3.15)
where ∇ 2 = ∇ • ∇ = ∆ is the Laplacian operator. Note that this energy can also been seen as an approximation of:

D(z [j-.5] ) R T ∇ 2 p -RT ∇ 2 p 2 K b dxdy, ( 3
E vol (z [j-.5] ) = D(z [j-.5] ) k v (|∇p| -|∇p|) 2 dx dy (3.17)
Note that the determinant is rotational invariant, hence it is not necessary to rotate the Jacobians.

As explained in Appendix A.1 integrating the volume energy (3.17) leads to:

E vol (z [j-.5] ) = πr 2 k v s 2 ∇ z c -s2 w wT ∇ z c 2 2 (3.18) + 1 2 πr 4 k v s 3 Ω u -s3 Ωu 2 (3.19) + 1 2 πr 4 k v s 3 Ω v -s3 Ωv 2 .
(3.20) Figure 3.7: Being built on VIPER primitives, our physics engine can simulate soft-body deformation and dynamic interactions between hundreds of models in real-time. A peculiarity of our engine is that both collision and simulation are executed on the same geometry; see our video in the additional material.

Optimization

To approach the non-linear optimization problem in (3.9), we linearize the inertia and non-linear elastic terms, and then solve the optimization iteratively in a Gauss-Newton fashion. To warm start the optimization, we first compute a prediction step by ignoring the elastic potentials and by solely minimizing the inertia term -this simply provides an initial guess. We then compute a (set of) correction steps that also include the elastic potentials.

Prediction step. By denoting by θ the angles parametrizing the rotation matrix and I is the second moment of area of a disc in world-space, the predictions for the different degrees of freedom of our model are computed as:

ĉ = c + h ċ + h 2 πr 2 m ξ ext , (center prediction) (3.21) θ = θ + h θ + I -1 h 2 s 2 m τ ext , (frame prediction) (3.22) ŝ = s + h ṡ + 2h 2 πr 4 m γ ext , (scale prediction) (3.23)
where ξ ext is the sum of the external forces which act on the disc, τ ext is the sum of the external torques, and γ ext is a quantity which can be seen as the counterpart of the total external torque measuring the sum of the external forces projected on the position vectors; see Appendix A.2 for more details. Note that the center and frame predictions are similar to the rigid-body equations of motion for a stretched disc. On top of these equations, we get a scale prediction describing how the scale of the disc is affected by the velocity and the external forces.

Correction steps. We then compute a set of correction steps including both inertial/elastic terms. We define

X = [c T [0] , s [0] , θ T [.5] , c T [1] , s [1] , θ T [1.5] , • • • ] T (3.24)
as the vector containing all the degrees of freedom, and λ as the vector of Lagrange multipliers.

K is a block diagonal matrix containing the stiffness parameters multiplied by the length of the piecewise elements, A is a block diagonal matrix stacking the inertia weights multiplied by the length of the piecewise elements, and

W(X) = [W 1 (X), W 2 (X), • • • ] T (3.25)
stacks the potential energy functions. Denoting the iteration number with k, the the state is then updated as X k = X k-1 + ∆X and λ k = λ k-1 + ∆λ, where, as derived in Appendix A.3:

∆X = -h 2 A -1 ∇W(X k-1 ) T ∆λ, (3.26) ∆λ = h 2 ∇W(X k-1 ) T 2 A -1 +K -1 -1 W(X k-1 ) -K -1 λ k-1 . (3.27)
For realtime performance we opt for using an iterative linear system solver such as block Jacobi or Gauss-Seidel. The update for the i-th constraint is:

∆λ i = β h 2 W i (X) T 2 A -1 + K -1 i -1 W i (X) -K -1 i λ i (3.28)
where we dropped the superscripts to improve readability, and β is a relaxation parameter. Note that A being a block diagonal matrix with block of size at most 3 × 3, A -1 can be efficiently computed. Note that Equation 3.28 is a generalization of the XPBD update [START_REF] Macklin | XPBD: Position-based Simulation of Compliant Constrained Dynamics[END_REF] derived for our volume preserving rod model.

Real-time physics engine

We implemented our real-time solver on a GPU by leveraging the Thrust framework provided by the CUDA library, which we execute on a single NVIDIA GTX 1080 graphics card. In our engine, any volumetric object is modeled as a collection of tapered capsule primitives, or "pills" (for brevity), while the floor is represented as a simple halfplane constraint. During simulation, at each time step, we start by first animating kinematic objects (e.g. bones). We then compute the inertial predictions, and perform collision detection (Sec. 3.7.1) to generate collision resolution constraints (Sec. 3.7.2). We then solve for all constraints using a Jacobi solver: constraints are computed in parallel, and the resulting positional displacements are averaged out by a reduction.

The transformations of VIPERs can then be skinned to any surface mesh model (Sec. 3.7.4).

Our model provides a viable alternative to cages as a compact control structure for soft-body deformations. We demonstrate this in Figure 3.7, where we rigged a simple octopus character using rods, and solve for collisions and soft-body deformations in real time. As shown in the accompanying video, our non-optimized prototype achieves real-time performance (≈ 7ms sim, ≈ 6ms render) for scenes containing up to 100 octopuses, each rigged with 37 pills, whose deformation is skinned to a triangular surface mesh of ≈ 13k faces. Note how for robust collision detection we need far fewer pills than the volumetric particles used in [START_REF] Macklin | Unified particle physics for real-time applications[END_REF].

Collision detection

To detect collisions between physical primitives, we adopt the approach presented by [START_REF] Green | Particle simulation using CUDA[END_REF] popularized in the context of GPU particle fluid simulation. Towards this goal, we approximate each pill by its bounding sphere, and (conservatively) detect collisions to be later handled in the resolution loop detailed in Section 3.7.2. Collisions are detected with the assistance of a uniform grid with cell width chosen as the diameter of the largest bounding sphere, such that collisions between spheres centered in non-neighboring cells are impossible. Extending [START_REF] Green | Particle simulation using CUDA[END_REF], we also count the spheres in the neighboring cells of each particle, and use this information to construct in parallel a list of all potential collisions. This is in contrast to the original algorithm which for each sphere processes neighbors in series, and therefore degrades to a partially serial algorithm in cases where many particles fall into a single grid cell. Further details regarding this process are provided together with an executable 1D example in the form of a Jupyter notebook in our additional material.

Collision handling

In a generic physics engine one would implement collision detection/response between any pair of available proxies. For the sake of efficiency, in our framework we only tackle two collision proxies:

(kinematic) half-planes and (dynamic) pills. Within the combinatorial set of collision pairs, the main challenge is pill-to-pill collisions. In what follows, we first compute the meta-parameters of the collisions, that are then resolved in the optimization via PBD constraints [START_REF] Müller | Position based dynamics[END_REF].

Collision metadata. Given two pills P a and P b , each modeled as a pair of spheres, for example,

P a = {(c a 1 , r a 1 ), (c a 2 , r a 2 )
}, the fundamental queries we need to answer are: 1 is there a collision? 2 what is the collision point/normal? 3 what is the inter-penetration amount? As typical in efficient collision resolution, we introduce a single PBD constraint modeling collision forces corresponding to the largest pill-to-pill inter-penetration. In our solution, we leverage the geometric structure of the problem: 1 a pill can be interpreted as the union of infinitely many spheres whose position and radii are linearly interpolated between its endpoints; 2 the largest inter-penetration corresponds to the inter-penetration between any pair of spheres, one in pill P a and one in pill P b . By first defining the LERP operator as L(x 1 , x 2 , γ) ≡ (1 -γ) x 1 + γ x 2 , the interpolated sphere (c(γ), r(γ)) is derived by LERP'ing the endpoint quantities as c(γ) = L(c 1 , c 2 , γ) and r(γ) = L(r 1 , r 2 , γ). The largest inter-penetration is then given by the solution of the bivariate optimization problem:

arg min α,β c a (α) -c b (β) 2 -(r a (α) + r b (β)) (3.29)
Because a closed-form operator Π b (α) providing the barycentric coordinate of the closest-point projection of a point c a (α) onto the pill P b is available (see Appendix A.4) we can further simplify this problem into a scalar optimization problem:

arg min α c a (α) -c b (Π b (α)) 2 -(r a (α) + r b (Π b (α))) (3.30)
which we solve by Dichotomous Search [START_REF] Antoniou | Practical Optimization: Algorithms and Engineering Applications[END_REF] with a fixed number of iterations (set to 10 in our engine).

Collision constraints.

Having detected a collision between two pills P a and P b , and having computed α * (and hence β * = Π b (α * )) by solving (3.30), we can express a constraint that correlates radii and positions on the two pills in order to resolve the collision in a least squares sense:

E collision = ( c a (α * ) -c b (β * ) 2 -r a (α * ) -r b (β * )) 2 .
(3.31)

Scale-invariant shape matching -"Bundling"

To represent more complex geometry than individual rods, such as that in Figure 3.2, we can gather a collection of rods in a cross-section, and introduce a constraint to explain their joint deformation. We employ the assumption that muscle fibers in a muscle cross-section contract isotropically. Indexing the rods in a cross-section by i, our rod deformation model can be expressed as a similarity transform

T i =   s i R i c i 0 1   . (3.32)
As illustrated in Figure 3.9, for each muscle cross-section we define a scale invariant shapematching energy measuring the deviation of the current rod deformations T i from a global similarity transform T * of the rest deformations Ti as

E shape = i T * Ti -T i 2 2 . (3.33)
We treat this energy as a hard constraint by finding the optimal T * and setting T i = T * Ti as a post-processing step after few iterations. The optimal T * can be computed following the derivation in [START_REF] Umeyama | Least-squares estimation of transformation parameters between two point patterns[END_REF]. The optimal rotation R * can be found by solving

R * = arg min R∈SO(3) R -Σ i [s i R i c i -µ][s i Ri ci -μ] T 2 2 , ( 3.34)
rest configuration (and constraints) traditional shape-matching scale-invariant shape-matching Figure 3.9: Scale-invariant shape matching -Shape matching can recover a rigidly transformed configuration, while our model allows for a null-space that includes uniform scale. We employ this model as we work under the assumption that muscle fibers in a muscle cross-section contract isotropically.

muscle at rest (not activated) traditional shape-matching scale-invariant shape-matching where µ = 1 n Σ i c i . We compute the optimal rotation R * using the robust approach presented in [START_REF] Müller | A robust method to extract the rotational part of deformations[END_REF]. The optimal scale can be computed as

s * = i sum(R * [s i Ri ci -μ] • [s i R i c i -µ]) i sum([s i Ri ci -μ] • [s i Ri ci -μ])
, (3.35) where sum(•) adds all entries of the matrix and • is the Hadamard product. Finally, the optimal translation can be derived as

c * = µ -s * R * μ.
(3.36)

Skinning VIPERs to surface deformation

Our VIPER rods can also be employed as a volumetric proxy driving the deformation of a surface mesh; see Figure 3.1 and Figure 3.8. In particular, we blend the relative transformations between deformed and rest pose configurations via linear blend skinning. To achieve this, we utilize a modified version of the LBS weight computation by [START_REF] Thiery | Sphere-meshes: shape approximation using spherical quadric error metrics[END_REF]. In that paper, weights were computed by fairing an initial assignment where each vertex was fully attached to the nearest element. Instead, we modify this assignment to begin with weights for each vertex that are proportional to the inverse square-distance from the vertex to the surface of each pill. This resolves cases where multiple surfaces are at equal distance, and the nearest neighbor is multiply-defined.

Energy implementation

The energies derived in Section 3.5 have been derived using continuous operators. To implement these energies we approximate ∇ z c(z [j+.5] ), ∇ z s(z [j+.5] ), and ∇ 2 z s(z [j+.5] ) using finite difference such that

∇ z c(z [j+.5] ) = l -1 [j] c(z [j+1] ) -c(z [j] ) , (3.37)
∇ z s(z [j+.5] ) = l -1 [j] s(z [j+1] ) -s(z [j] ) , (3.38) ∇ 2 z s(z [j] ) = l -1 [j] s(z [j+1] ) -s(z [j] ) (3.39) -l -1 [j-1] s(z [j] ) -s(z [j-1] ) .
For simplicity of the derivations we used angles θ to parametrize the rotation R . However, our implementation uses quaternions. For small angles the corresponding quaternion is

Q = [ θ 2 T , 1] T .
We also use quaternions to represent rotations and approximate the Darboux vector using

Ω(z [j] ) = 4(l [j-1] + l [j] ) -1 Im( Q[j-.5] Q [j+.5] ), (3.40)
where Im(•) gives the imaginary part of a quaternion.

Anatomical modeling and simulation

We now describe how our rods can be used to model complex anatomical structures such as bones and muscles; see Figure 3.1. For the former, we use a simplified version of [START_REF] Thiery | Sphere-meshes: shape approximation using spherical quadric error metrics[END_REF] where only pill primitives are used -this primarily allows us to reduce the complexity of the collision detection/resolution codebase. We then detail the conversion of digital models of muscles into VIPERs in Section 3.8.1, and describe a few nuances about the simulation of their motion in Section 3.8.2.

Muscles to rods conversion -"Viperization"

We developed a (weakly-assisted) technique to convert a traditional muscle model into a collection of K rods. As outlined in Figure 3.11, our process involves several phases. We begin by computing a volumetric discretization of the muscle's surface. We then ask the artist to paint annotations on the surface marking the start (sources) and the end (sinks) of the muscle. A harmonic solve alike the one describe in Choi and Blemker [2013] is the executed to compute a field that smoothly varies in the [0, 1] range along the muscle's length. We then sample M iso-levels of this field to be surfaces containing the M vertices of each of the K generated VIPERs. Within each iso-level we require: 1 VIPERs to have the same radii, and 2 to be distributed uniformly on the slice (iso-surface). To obtain this, we perform a restricted centroidal Voronoi diagram of K points on each slice [START_REF] Botsch | Projective dynamics: Fusing constraint projections for fast simulation[END_REF], while simultaneously penalizing the length of each rod. We alternate this variational optimization with a discrete one that re-assigns spheres to different rods in order to minimize the sum of rod lengths. This process, which we refer to as "combing", starts from one end of the muscle, and executes in order M -1 instances of minimum-cost bipartite matching (which we solve in polynomial time via the Hungarian algorithm), where the pairwise costs are the K 2 euclidean distances across rod nodes in two adjacent slices. A few examples of the results of this process are illustrated in Figure 3.1 and Figure 3.12.

Muscle simulation

Once a muscle is viperized as described in Section 3.8.1, rod centers within the same cross-section are connected via the bundling constraints in Equation 3.32. Every rod also obeys the constraints described in Section 3.5. The endpoints of rods are kinematically attached to bones, and intramuscle collisions are disabled, while inter-muscle collisions are detected and resolved as described respectively in Section 3.7.1 and Section 3.7.2. Muscles can also be activated (fiber contractions generating a stronger force, producing a change of shape at constant muscle length and volume) by inserting internal forces, or even simply shortening the length of fibers resulting in the bulging effects illustrated in Figure 3.10. We also speed-up the solver convergence during fast motion by exploiting the availability of bone transformations. In more detail, each muscle particle has two skinning weights corresponding to the two bones the muscle is attached to. At the beginning of each frame we use the transform of each bone relative to their last frame's transforms to initialize the displacement of the particle using LBS, and later refine this via simulation. The mesh-to-VIPER conversion process. Given source/sink constraints, we compute a harmonic function in the volume, and extract a few discrete iso-levels. Within each of these, we execute a restricted CVD to place 5 elements of the same radii on these surfaces. We then execute a combinatorial optimization that connects samples across layers to produce minimal length curves.

Conclusions & future work

In this paper, we introduced a novel formulation of cosserat rods that considers local volume, and optimizes for its local conservation. The resulting position-based simulation is highly efficient, and is straightforward to implement on graphics hardware. We demonstrated how rod-bundling is a powerful representation for the modeling of volumetric deformation -and in particular for skeletal muscles. Rather than requiring artists to model from scratch, we also introduced an algorithm to procedurally generate VIPERs with minimal user interaction. Finally, by coupling the rod simulation to a surface mesh via skinning, our model can be thought of as a direct alternative to tetrahedral meshes and cages for real-time non-rigid deformation. Most importantly, our generalized rods formulation opens up a number of venues for future work, which we classify in three broad areas, as elaborated below.

Model generalization. While in our rods we discretized the skeletal curve with piecewise linear elements, it would be interesting to investigate whether the use of continuous curve models such as splines would be tractable -from both mathematical, as well as implementation standpoints.

Similarly to [Müller and Chentanez, 2011;Müller and Chentanez, 2011], our model could also be extended to model anisotropic volume deformations, that is, both the rest pose and deformed rod could have a non-circular cross section. Further, while in this paper we treated the modeling of non-circular cross-sections via bundling, the theory of medial axis [START_REF] Tagliasacchi | 3d skeletons: A state-of-the-art report[END_REF] tells us how any geometric object can be approximated via primitives formed as the convex-hulls of three-spheres -what [START_REF] Tkach | Sphere-meshes for real-time hand modeling and tracking[END_REF] called "wedges". Extending our volume-invariant rods models to volume-invariant wedges would provide an elegant generalization of our modeling paradigm.

Anatomical modeling. As highlighted by our supplementary video, rod primitives can be exploited for the efficient approximate modeling and simulation of complex structures. Nonetheless, the dynamics of the human body are the result of the complex interplay between muscle, fat, and the consequent deformation of skin. Enriching our model to also account for these factors would be an interesting extension. For example, rather than driving the muscle surface via skinning as in Section 3.7.4, one could represent a muscle as a controllable implicit blend [START_REF] Angles | Sketch-based implicit blending[END_REF], approximate fat as an elastic offset between muscles and skin, and simulate skin as an elastic surface whose vertices lie on a (potentially) user-controlled iso-level of the implicit function.

Further, while artistical editing of physically driven anatomical systems can be a difficult due to the complexity of simulation, our framework could immediately enable interactive modeling, in a similar fashion to what ZBrush/ZSphere currently provides for authoring static geometry. By extending the works in [START_REF] Tkach | Sphere-meshes for real-time hand modeling and tracking[END_REF][START_REF] Tkach | Online generative model personalization for hand tracking[END_REF], an efficient anatomical model for a particular user could also be constructed by fitting to RGBD data.

Optimization. Our solver has not yet directly leveraged the straightforward multi-resolution structure of rod geometry. More specifically, the curve parameterization of rods offers a domain over which designing prolongation/restriction operators needed for a geometric multi-grid implementation becomes straightforward. An orthogonal dimension for optimization would be to consider the existence of multi-resolution structures within the cross-sectional domain; see Figure 3.2. This could both be exploited in offering multi-scale interaction for artists in editing our rod models, as well as to produce level-of-details models for efficient simulation at scale. Finally, while we employed out-of-the-box geometry processing tools to convert a triangular surface mesh into a rod model, we believe fitting a fiber-bundle model to a given solid could be achieved without the (often finicky) conversion to volume/tetrahedral mesh, but rather as a direct optimization over fiber placements. Patches are then extracted via a sampler S, and then fed to a task-specific network T τ . The results are then fed to a task-specific loss L. In this example, the specific task is to re-synthesize the image as a super-position of (unknown, local) basis functions.

of sequential detection [START_REF] Eslami | Attend, Infer, Repeat: Fast Scene Understating with Generative Models[END_REF] or channel-wise approaches [Zhang et al., 2018c] are problematic to apply when the number of instances of the same object is large. Applications. To demonstrate the capabilities of MISTs, we evaluate our network on a variety of weakly-supervised multi-instance problems. Note how in some of these applications, the value of K is the only supervision signal we provide. We consider 1 the problem of recovering the basis functions that created a given texture, 2 the classification of numbers in cluttered scenes where the only supervision is the occurrence of these numbers.

Introducing

In summary, in this paper:

• we introduce the MIST framework for weakly-supervised multi-instance visual learning;

• we introduce a training method that allows the use of top-K approaches for end-to-end trainable architectures;

• we show that our framework can reconstruct images as parts, as well as detect/classify instances without any location supervision.

Related works

Attention models and the use of localized information have been actively investigated in the literature. Some examples include discriminative tasks such as fine-grained classification [START_REF] Sun | Multi-Attention Multi-Class Constraint for Fine-grained Image Recognition[END_REF], and pedestrian detection [Zhang et al., 2018a], and generative ones such as image synthesis from natural language [START_REF] Johnson | Image Generation from Scene Graphs[END_REF]. We now discuss a selection of representative works, and classify them according to how they deal with multiple instances.

Grid-based methods. Since the introduction of Region Proposal Networks (RPN) [START_REF] Ren | Faster R-CNN: Towards Real-Time Object Detection with Region Proposal Networks[END_REF], grid-based strategies have been used for dense image captioning [START_REF] Johnson | Densecap: Fully Convolutional Localization Networks for Dense Captioning[END_REF],

instance segmentation [START_REF] He | Mask R-CNN[END_REF], keypoint detection [START_REF] Georgakis | End-to-end Learning of Keypoint Detector and Descriptor for Pose Invariant 3D Matching[END_REF], multi-instance object detection [START_REF] Redmon | YOLO 9000: Better, Faster, Stronger[END_REF]. Recent improvements to RPNs attempt to learn the concept of a generic object covering multiple classes [START_REF] Singh | at 30fps: Decoupling Detection and Classification[END_REF], and to model multi-scale information [START_REF] Chao | Rethinking the Faster R-CNN Architecture for Temporal Action Localization[END_REF]. The multiple transformation corresponding to separate instances can also be densely regressed via Instance Spatial Transformers [START_REF] Wang | Geometry-Aware Scene Text Detection with Instance Transformation Network[END_REF], which removes the need to identify discrete instance early in the network. Unfortunately, all these method are fully supervised, as they require both class labels and object locations for training.

Heatmap-based methods. Heatmap-based methods have recently gained interest to detect features [START_REF] Yi | LIFT: Learned Invariant Feature Transform[END_REF][START_REF] Ono | Lf-Net: Learning Local Features from Images[END_REF][START_REF] Detone | Superpoint: Self-Supervised Interest Point Detection and Description[END_REF], find landmarks [Zhang et al., 2018c;[START_REF] Merget | Robust Facial Landmark Detection via a Fully-Conlolutional Local-Global Context Network[END_REF], and regress human body keypoint [START_REF] Tekin | Learning to Fuse 2D and 3D Image Cues for Monocular Body Pose Estimation[END_REF][START_REF] Newell | Stacked Hourglass Networks for Human Pose Estimation[END_REF].

While it is possible to output one heatmaps per type of point [Zhang et al., 2018c;[START_REF] Tekin | Learning to Fuse 2D and 3D Image Cues for Monocular Body Pose Estimation[END_REF], this still restricts the number of instances to one. [START_REF] Yi | LIFT: Learned Invariant Feature Transform[END_REF] re-formulates the problem based on each instance, but in doing so it introduces a non-ideal difference between training and testing regimes. Grids can also be used in combination to heatmaps [START_REF] Detone | Superpoint: Self-Supervised Interest Point Detection and Description[END_REF], but this results in an unrealistic underlying assumption of uniformly distributed detections in the image. Overall, heatmap-based methods excel when the "final" task of the network is generate a heatmap [START_REF] Merget | Robust Facial Landmark Detection via a Fully-Conlolutional Local-Global Context Network[END_REF], but are problematic to use as an intermediate layer in the presence of multiple instances.

Sequential inference methods.

Another way to approach multi-instance problems is to attend to one instance at a time in a sequential way. [START_REF] Gregor | DRAW: A Recurrent Neural Network For Image Generation[END_REF] proposes a recurrent network that processes only a small area at a time for both discriminative and generative tasks.

These sequential models have then been extended to localize and recognize MNIST digits in a cluttered image [START_REF] Ba | Multiple Object Recognition With Visual Attention[END_REF][START_REF] Eslami | Attend, Infer, Repeat: Fast Scene Understating with Generative Models[END_REF]. Overall, RNNs often struggle to generalize to sequences longer than the ones encountered during training, and while recent results on inductive reasoning are promising [START_REF] Gupta | Inductive Visual Localization: Factorised Training for Superior Generalization[END_REF], their performance does not scale well when the number of instances is large.

Knowledge transfer. To overcome the acquisition cost of labelled training data, one can transfer knowledge from labeled to unlabeled dataset. For example, Inoue et al. [START_REF] Inoue | Cross-Domain Weakly-Supervised Object Detection through Progressive Domain Adaptation[END_REF] train on a single instance dataset, and then attempt to generalize to multi-instance domains, while Uijlings et al. [START_REF] Uijlings | Revisiting Knowledge Transfer for Training Object Class Detectors[END_REF] attempts to also transfer a multi-class proposal generator to the new domain. While knowledge transfer can be effective, it is highly desirable to devise unsupervised methods such as ours that do not depend on an additional dataset.

Weakly supervised methods. To further reduce the labeling effort, weakly supervised methods have also been proposed. [START_REF] Wan | Min-Entropy Latent Model for Weakly Supervised Object Detection[END_REF] learns how to detect multiple instances of a single object via region proposals and ROI pooling, while Tang et al. [Tang et al., 2018b] proposes to use a hierarchical setup to refine their estimates. [START_REF] Gao | C-WSL: Coung-guided Weakly Supervised Localization[END_REF] provides an additional supervision by specifying the number of instances in each class, while Zhang et al.

[ Zhang et al., 2018b] localizes objects by looking at the network activation maps [START_REF] Zhou | Learning Deep Features for Discriminative Localization[END_REF][START_REF] Selvaraju | Grad-CAM: Visual Explanations from Deep Networks via Gradient-based Localization[END_REF]. However, all these method still rely on region proposals from an existing method, or define them via a hand-tuned process.

MIST Framework

A prototypical MIST architecture is composed of two trainable components: 1 the first module receives an image as input and extracts a collection of patches, at image locations and scales that are computed by a trainable heatmap network H η with weights η; see Section 4.4. 2 the second module processes each extracted patch with a task-specific network T τ whose weights τ are shared across patches, and further manipulates these signals to express a task-specific loss L task ; see Section 4.5. The two modules are connected through non-maximum suppression on the scale-space heatmap output of H η , followed by a top-K selection process to extract the parameters defining the patches, which we denote as E K . We then sample patches at these locations through bilinear sampling S and feed them the second module.

The defining characteristic of MIST architectures is that they are quasi-unsupervised: the only strictly required supervision is the number K of patches to extract. The training of MIST architectures is summarized by the optimization:

arg min τ,η L task (T τ (S(E K (H η (I))))) (4.1)
where τ, η are the network trainable parameters. Note how in the expression above E K is nondifferentiable, thus making (4.1) unapproachable by back-propagation for training.

Example task. In Figure 4.1, we illustrate an example of a MIST architecture for image reconstruction. In more details, the task is to understand how to re-synthesize the image as a super-position of K spatially localized basis functions. Note that, while this task is quasiunsupervised, it presents several joint challenges as it needs to estimate: 1 an unknown shared low-dimensional set of latent bases; 2 where to place instances from this latent space;

Patch extraction

We extracts a set of K (square) patches that correspond to "important" locations in the image -where importance is a direct consequence of L task . The localization of such patches can be computed by regressing a 2D heatmap whose top-K peaks correspond to the patch centers.

However, as we do not assume these patches to be equal in size, we regress to a collection of heatmaps at different scales. To limit the number of necessary scales, we use a discrete and sparse scale-space, while resolving for intermediate scales by weighted interpolation.

Multiscale heatmap network -H η . Our multiscale heatmap network is inspired by LF-Net [START_REF] Ono | Lf-Net: Learning Local Features from Images[END_REF]. We employ a fully convolutional network with (shared) weights η at multiple scales, indexed by s = 1 . . . S, on the input image I. The weights η across scales are shared so that the network cannot implicitly favor a particular scale. To do so, we first downsample the image to each scale I s , execute the network H η on it, and finally upsample to the original resolution. This process generates a multiscale heatmap tensor h = {h s } of size H × W × S where h s = H η (I s ), and H is the height of the image and W is the width. For the convolutional network we use 4 ResNet blocks [START_REF] He | Delving Deep into Rectifiers: Surpassing Human-Level Performance on Imagenet Classification[END_REF], where each block is composed of two 3 × 3 convolutions with 32 channels and relu activations without any downsampling. We then perform a local spatial softmax operator [START_REF] Ono | Lf-Net: Learning Local Features from Images[END_REF] with spatial extent of 15 × 15 to sharpen the responses.

Extracting patch location and scale. We first normalize the heatmap tensor so that it has a unit sum along the scale dimension:

h s = H η (I s ) s H η (I s ) + ǫ , ( 4.2)
where ǫ = 10 -6 is added to prevent divisions by zero. We then compute the top-K spatial locations across all scales, to obtain the image-space coordinates of patch centers:

{(x k , y k )} = TOPK(max s h s ). (4.3)
Note that a direct extraction of K maxima is possible because the aforementioned local spatial softmax performs a localized non-maximal suppression. The corresponding scale is computed by weighted first order moments [START_REF] Suwajanakorn | Discovery of Latent 3D Keypoints via End-To-End Geometric Reasoning[END_REF], where the weights are the responses in the corresponding heatmaps: Note also that our extraction process uses a single heatmap for all instances that we extract.

s k = s h s (x k , y k ) weights s. ( 4
By contrast, existing heatmap-based methods [START_REF] Eslami | Attend, Infer, Repeat: Fast Scene Understating with Generative Models[END_REF]Zhang et al., 2018c] typically rely on heatmaps dedicated to each instance, which is problematic when an image contains two instances of the same class. Conversely, we restrict the role of the heatmap network H η to find the "important" areas in a given image, without having to distinguishing between classes, hence simplifying learning.

Patch sampling. As a patch is uniquely parameterized its location and scale

x k = (x k , y k , s k ),
we can then proceed to sample its corresponding tensor via bilinear interpolation [START_REF] Jaderberg | Spatial Transformer Networks[END_REF]:

{P k } = S (I, {x k }) . ( 4

.5)

Comparison to LF-Net. Note that differently from LF-Net [START_REF] Ono | Lf-Net: Learning Local Features from Images[END_REF], we do not perform a softmax along the scale dimension. The scale-wise softmax in LF-Net is problematic as the computation for a softmax function relies on the input to the softmax being unbounded. For example, in order for the softmax function to behave as a max function, due to exponentiation, it is necessary that one of the input value reaches infinity (i.e. the value that will correspond to the max), or that all other values to reach negative infinity. However, at the network stage where softmax is applied in [START_REF] Ono | Lf-Net: Learning Local Features from Images[END_REF], the score range from zero to one, effectively making the softmax behave similarly to averaging. Our formulation does not suffer from this drawback.

Task-specific networks

We now introduce two instances of the MIST architectures corresponding to different applications.

We keep the heatmap network and the extractor architectures the same, and only change the taskspecific network, as well as the loss used for supervision. In particular, we consider a reconstruction problem in Section 4.5.1, and a classification problem in Section 4.5.2. Further implementation details can be found in Section 4.8.

Image reconstruction

As illustrated in Fig. 4.1, for image reconstruction we append our patch extraction network with a shared auto-encoder for each extracted patch. We can then train this network to reconstruct the original image by inverting the patch extraction process, and forming the task specific loss to be the ℓ 2 norm between the input image and the reconstructed image. Specifically, we introduce the inverse sampling operation S -1 (P i , x i ), which starts with an image of all zeros, and places the patch P i at x i . We then add all the images together to obtain the reconstructed image, overall expressing the following task loss: h

L task = I - i S -1 (P i , x i ) 2 2 . ( 4
← E -1 K ({x * k }) 8: η ← Optimize H η with L lift 9:
end for 10: end function Overall, the network is designed to jointly model and localize repeating structures in the input signal. Regressing shared basis functions can be related to non-local mean processes [START_REF] Buades | A Non-Local Algorithm for Image Denoising[END_REF], as a model for the input signal is created by agglomerating the information in scattered spatial instances. Our task architecture is also related to transforming auto-encoders [START_REF] Hinton | Transforming Auto-Encoders[END_REF], where the difference is that in this previous work a single instance is present in the image, and that they provide the ground truth transformations as a supervision.

Multiple instances classification

By appending the patch extraction module with a classification network we can realize an architecture for multiple instance learning. For each extracted patch P k we apply a shared classifier network to output ŷk ∈ R C , where C is the number of classes. In turn, these are then converted into probability estimates by the transformation pk = softmax(ŷ k ). By denoting the one-hot ground-truth labels of instance l in the image as y l , we define the multi-instance classification loss as

L task = H 1 L L l=1 y l , 1 K K k=1 pk , ( 4.7)
where H (•) denotes cross entropy and L is the number of instances in the image. Note here that we do not provide supervision about where each class instances are, yet the detector network will automatically learn how to localize the content with minimal supervision.

Training MISTs

While it is technically possible to back-propagate through MIST architectures, the gradients would only flow through the spatial regions corresponding to the K selected keypoints -this results in a training process that ignores locations away from the selection. This is a fundamental issue, as, in Differentiable top-K via lifting. The introduction of auxiliary variables (i.e. lifting) to simplify the structure of an optimization problem has proven effective in a range of domains ranging from registration via ICP [START_REF] Taylor | Efficient and precise interactive hand tracking through joint, continuous optimization of pose and correspondences[END_REF], to efficient deformation models [START_REF] Sorkine | As-rigid-as-possible surface modeling[END_REF], and robust optimization [START_REF] Zach | Descending, Lifting or Smoothing: Secrets of Robust Cost Opimization[END_REF]. To simplify our training optimization, we start by decoupling the heatmap tensor from the optimization (4.1) by introducing the corresponding auxiliary variables h, as well as the patch parameterization variables {x k } that are extracted by the top-K extractor:

arg min η,τ, h,{x k } L task (T τ (S({x k }))) (4.8) s.t. {x k } = E K ( h) (4.9) h = H η (I) (4.10)
We then relax (4.10) to a least-squares penalty:

arg min η,τ, h,{x k } L task (T τ (S({x k }))) + h -H η (I) 2 2 (4.11) s.t. {x k } = E K ( h) (4.12)
and finally approach it by alternating optimization:

arg min τ,{x k } L task (T τ (S({x k }))) (4.13) arg min η h -H η (I) 2 2 (4.14)
where h has been dropped as it is not a free parameter: it can be computed as h = E -1 K ({x k }) after the {x k } have been optimized by (4.13), and as h = H η (I) after η have been optimized by (4.14).

To accelerate training, we further split (4.13) into two stages, and alternate between optimizing τ and {x k }. In particular, multiple optimization iterations of {x k } are executed to allow keypoints to displace faster during training. The summary for the three stage optimization procedure is outlined in Alg. 1: 1 we optimize the parameters τ with the loss L task ; 2 we then fix τ , and refine the positions of the patches {x k } for T iterations with L task . 3 with the optimized patch positions {x * k }, we invert the top-K operation by creating a target heatmap h, and optimize the parameters η of our heatmap network H using ℓ 2 distance between the two heatmaps, L lift = h -H η (I) 2 2 . Notice that we are not introducing any additional supervision signal that is tangent to the given task.

Generating the target heatmap -E -1 K ({x k }). For creating the target heatmap h, we create a tensor that has zeros everywhere except for the positions corresponding to the optimized positions.

However, as the optimized patch parameters are no longer integer values, we need to quantize them with care. For the spatial locations we simply round to the nearest pixel, which at most creates a quantization error of half a pixel, which does not cause problems in practice. For scale however, simple nearest-neighbor assignment causes too much quantization error as our scale-space is sparsely sampled. We therefore assign values to the two nearest neighboring scales in a way that the center of mass would be the optimized scale value. That is, we create a heatmap tensor (bottom) the hard multi-instance setup. We compare the output of MISTs to grid, channel-wise, Eslami et al. [START_REF] Eslami | Attend, Infer, Repeat: Fast Scene Understating with Generative Models[END_REF] and Zhang et al. [Zhang et al., 2018c].

that would result in the optimized patch locations when used in forward inference.

Results and evaluation

To demonstrate the effectiveness of our framework we evaluate two different tasks. We first perform a quasi-unsupervised image reconstruction task, where only the total number of instances in the scene is provided. We then show that our method can also be applied to weakly supervised multi-instance classification, where only image-level supervision is provided. Note that, unlike region proposal based methods, our localization network only relies on cues from the classifier, and both networks are trained from scratch.

Image reconstruction

From the MNIST dataset, we derive two different scenarios. In the "MNIST easy" dataset, we consider a simple setup where the sorted digits are confined to a perturbed grid layout; see 4.2 (top). Specifically, we perturb the digits with a Gaussian noise centered at each grid center, with a standard deviation that is equal to one-eighths of the grid width/height. In the "MNIST hard" dataset, the positions are randomized through a Poisson distribution [START_REF] Bridson | Fast Poisson Disk Sampling in Arbitrary Dimensions[END_REF],

as is the identity, and cardinality of each digit. Note how we allow multiple instances of the same digit to appear in this variant. As expected, both these datasets contain a training and testing subsets, and the testing portion is never seen at training time.

Comparison baselines. We compare out method against four baselines: 1 grid we setup a grid of keypoints, and apply the same auto-encoder architecture as MIST to reconstruct the input image; 2 in the channel-wise variant we use the same heatmap network, except for the last convolutional layer giving K channels as output, where each channel is dedicated to an interest point. Their locations are obtained through a channel-wise soft argmax as in [Zhang et al., 2018c].

We also use the same architecture for the auto-encoder as MIST; 3 the method of Eslami et al.

[ [START_REF] Eslami | Attend, Infer, Repeat: Fast Scene Understating with Generative Models[END_REF] is a sequential generative model. To generate nine digits, it is required for the framework is tasked to regress the underlying impulse process, and reconstruct the Gabor kernels so that when the two are combined, we can reconstruct the original image. Gabor kernels, even though in the training images they are heavily overlapped. Further, note that the number of instances detected is significantly larger than that possible with other methods.

Multiple instance classification

To test our method in a multiple instance classification setup, we rely on the MNIST hard dataset.

We compare our method to channel-wise, as other baselines are designed for purely generative tasks.

To evaluate the accuracy of the models, we compute the intersection over union (IoU) between the ground-truth bounding box and the detection results, and assign it as a match if the IoU score is over 50%. We report the number of correctly classified matches in Table 4.2. Our method clearly outperforms the channel-wise strategy. A few qualitative results are illustrated in Figure 4.5. Note that even without direct supervision on the digits locations, our method correctly localizes them.

Conversely, the channel-wise strategy fails to learn. This is because multiple instances of the same digits are present in the image. For example, in the example Figure 4.5 (right), we have two number sizes, zeros, and nines. This prevents any of these digits from being detected/classified properly by a channel-wise approach.

Implementation details

Auto-encoder network. The input layer of the autoencoder is 32x32xC where C is the number of color channels. We use 5 up/down-sampling levels. Each level is made of 3 resnet blocks and each resnet block uses a number of channels that doubles after each downsampling step. Resnet blocks uses 3x3 convolutions of stride 1 with relu activation. For downsampling we use 2D max pooling with 2x2 stride and kernel. For upsampling we use 2D transposed convolutions with 2x2 stride and kernel. The output layer uses a sigmoid function. We use layer normalization before each convolution layer.

Classification network. We re-use the same architecture as encoder for first the task and append a dense layer to map the latent space to the score vector of our 10 digit classes.

Conclusion

In this paper, we introduced the MIST framework for multi-instance image reconstruction and classification. Both these tasks are based on localized analysis of the image, yet we train the network without providing any localization supervision. The network learns how to extract patches on its own, and these patches are then fed to a task-specific network to realize an end goal. While at first glance the MIST framework might appear non-differentiable, we show how via lifting they can be effectively trained in an end-to-end fashion. We demonstrated the effectiveness of MNIST by introducing a variant of the MIST dataset, and demonstrating compelling performance in both reconstruction and classification. We also show how the network can be trained to reverse engineer a procedural texture synthesis process.

Chapter 5

Conclusion

Summary

In Chapter 2, we have presented a method to reverse-engineer how primitives blend together. We have introduced a template for blending operators that can mimic any operator from previous work but also completely new ones. This template can be optimized to fit sample points in the operator domain. This enables the interactive design of operators that can be applied by users in composition modeling without necessitating any technical knowledge.

In Chapter 3, we have introduced a new volume-preserving primitive for physical simulations. In order to efficiently simulate it, we have extended the position-based dynamics simulation framework by adding an extra degree of freedom. To do so, we have carefully derived the equations of the updates and constraints. We have demonstrated that this primitive is an effective collision proxy.

Additionally, this primitive can be used as a volumetric proxy for surface deformation. Finally, we have shown that this simple primitive can be useful to simulate complex objects such as muscles.

In Chapter 4, we have introduced a novel unsupervised Deep Learning framework to learn and detect recurring primitives. At the core of our contribution is our multi-steps training process that is required to train our network despite the non-differentiable top-k operation. We have successfully applied our method on multi-instance reconstruction and classification tasks, outperforming the state-of-the-art unsupervised methods.

Future Work

Joint optimization of primitives and composition functions. In Chapter 2, a method to extract the function by which some primitives are combined was presented. This assumes that the shape of the primitives and their pose are known. But this information is rarely available. Most three-dimensional models are just specified in terms of their final shape, not by the primitives that compose them. Of course, there has been some work on extracting primitives but it would be interesting to jointly optimize for the primitives and the composition functions.

We aim at finding the unknowns x t = [ c T t st θ T t ] T that minimize (A.24) giving us the prediction update. We denote by θ the angles parametrizing the rotation matrix. Because the deformation function p t is non linear, i.e., due to the rotational degrees of freedom, we rely on a Gauss-Newton iterative scheme for the minimization. We linearize (A. where k is the iteration number, and ∆x t = [ ∆c T t ∆st ∆θ T t ] T . At each iteration we minimize (A.25), and then apply the update x k+1 t = x k t + ∆x k t , where we initialize x 0 t = x t-1 . The matrix A can be written as A = [ I3×3 R k t q -s k t [R k t q]× ], where [•] × is a cross product skew-symmetric matrix. The vector b is defined as b = p t-1 + h ṗt-1 + h 2 m f ext -p k t . To compute the prediction updates we will use a single iteration of Gauss-Newton so b = h ṗt-1 + h 2 m f ext as p 0 t = p t-1 . We will now drop the superscripts and the subscripts to improve readability. As Equation A.25 is quadratic we can be find its minimum by solving the normal equation

4. 1

 1 Figure 1.1: Examples of recurring primitives in nature. (a) A pine cone is composed on dozens of instances of the same primitive shape. (b) The skin of a snake exhibits a repetition of small diamond-shaped scales. (c) Bees build honeycomb structures that are made of a single volumetric primitive. (Sources: pratheep.com, Care SMC, Dominic Heard)

Figure 1 . 2 :

 12 Figure 1.2: Examples of fractals in nature. (a) The Romanesco Broccoli has a fractal-like shape with 4 levels of subdivisions. (b) The branches of a tree split in smaller versions of themselves. (c) A stellar dendrite snowflake that is composed of a central hexagon whose corners spawn another hexagon each. Each secondary hexagon spawns small hexagons at each corner. (Sources: Albrecht Fietz, Sunil Patel, SnowCrystals.com)

Figure 1 . 4 :

 14 Figure 1.4: Visualization of features maps of a convolutional neural network. (left) The first layers detect edges and intensity contrasts. (middle) The mid-level layers build abstraction upon the previous layers and respond to textures and patterns. (right) The high-level layers identify the present of objects they are trained for in the mid-level feature maps. (Source: [LeCun and Ranzato, 2013])

Figure 2

 2 Figure 2.1: A user's 2D sketches (bottom-left) are used to exemplify desired ways in which implicit functions should be composed together. From these, our algorithm automatically derives new custom gradient-based composition operators (bottom-right). These can then be applied to combine any 3D (or 2D) implicit model (top) replicating the user's intentions, and including effects such as contacts, bulging deformation, or smooth blends.

Figure 2 . 2 :

 22 Figure 2.2: Proper design of an implicit composition operator can achieve classical CSG operations such as union, difference, intersection, as well as their smooth variants -i.e. blending.

Figure 2 .Figure 2 . 3 :

 223 Figure 2.3: A composition operator defined as a function of solely the values of the two input objects presents undesirable bulge and tunnel artifacts.

Figure 2 . 3 .Figure 2 . 4 :

 2324 Figure 2.4: Through inclusion of the gradients angle θ, an operator that is capable of resolving the shortcomings of implicit blending can be designed.

Figure 2 . 5 :

 25 Figure2.5: Given a pair of input primitives (two overlapping circles), the user sketches the desired resulting shape (a contact surface and a bulging effect). (a) The system generates a set of samples in the operator space D. (b) An operator template is fitted to the samples. (c) A dense regular sampling of the operator g is generated. The resulting operator, if applied to the initial pair of primitives, produces the desired behavior and can be applied wherever this effect is desired.

Figure 2

 2 The field has a value which decreases with the distance from the skeleton, according to a given falloff function; see[START_REF] Marschner | Fundamentals of computer graphics[END_REF]. The area where the field function has values |f (x) -0.5| < 0.5 is denoted the support of the implicit object. Outside its support, the field function equals zero and the primitive has no influence on the composition operations. The support is compact, i.e. bounded; seeFigure 2.7.

Figure 2 . 6 :

 26 Figure 2.6: Operator design feedback loop: because our pipeline (Fig. 2.5) has interactive response times, our system allows progressive refinement of the operator though successive strokes. Undesired blending artifacts (right) are corrected until a final operator is constructed which is capable of reproducing the desired behavior.

Figure 2 .

 2 Figure 2.7: Top: a sample s is drawn in the intersection of the supports of the two primitives (bold colored lines) generated by the skeletons (black dots). Lower-left: a zoom-in around s: the gradients of the two field functions f a and f b form an angle θ. Lower-right, the corresponding sample s in the operator's domain D (see Eq. 2.5).

Figure 2

 2 Figure 2.8: Our template is a surface P is made of two bi-quadratic patches, depicted with orange and violet control points -red ones are shared, creating a seam at the junction. This template is able to represent a variety of useful composition operators; see Section 2.6.2.

Figure 2 .

 2 Figure 2.8c) and asymmetric contact (see Figure 2.8d) amongst many others; see Figure 2.21.

Figure 2 .

 2 Figure 2.21efg. Our template can easily reproduce this situation by making the two patches partially coincide, thus realizing a non-manifold operator.

  of influence, so f c should exactly reproduce the values of f b (p) to ensure C 0 continuity. Analogous considerations apply to the f b (p) = 0, f a (p) = 1 and f b (p) = 1 boundary planes of D, leading to the constraints:

Figure 2 .

 2 Figure 2.10: (a) Symmetric contact operators a-la[START_REF] Cani | An implicit formulation for precise contact modeling between flexible solids[END_REF] can be constructed through the enforcement of symmetry in D. We introduce asymmetric contact operators which can model phenomena as: (b) a rubber ball hitting a concrete wall, or (c) a steel ball hitting a sheet of softer metal.

Figure 2 .

 2 Figure 2.11: Our sketched implicit operators are easily controllable by the artist. In this figure, we demonstrate several variants of two types of operators: (top) the smooth-union from[START_REF] Gourmel | A gradient-based implicit blend[END_REF], as well as the (bottom) implicit-contact from[Vaillant et al., 2013a].

Figure 2 .

 2 Figure 2.12: We define two symmetric contact operators through the sketches in the first column, generating g a and g b resulting in the composition visualized in the second column. By leveraging the consistency in the operator's parameterization, other operators can be generated as g k = (1 -k)g a + kg b ; for k ∈ (0, 1) we obtain operators whose behavior is intuitively interpolated (top), while extrapolation can be obtained for values of k outside of this range (bottom).

Figure 2 .

 2 Figure 2.13: Modeling quasi-biological structures such as a procedural curly-tree (top) and an oak leaf (bottom). We show the input primitives and the artist sketches (a), as well as the composition result (b). The leaf design starts with the placement of a few skeletal branches, from which the leaf boundary (d) or a single sketch (f) can be used to design a corresponding blending (e,g).

Figure 2 .

 2 Figure 2.14: A rubber ball that elongates according to its velocity, and squishes when it comes into contact with the floor in an artist-defined fashion.

Figure 2 .

 2 Figure 2.16: The same sphere lying in proximity of a plane can produce contact deformations simulating both hydrophobic and hydrophilic deformations according to the type of asymmetric blending designed by the artist.

Figure 2 .

 2 Figure 2.18: (top) Multiple sketches on different primitives are used to synthesize a blending operator that can be used to produce font variations as well as an entire "infinite resolution" font family. (bottom) Sketching negative blends to apply a "serification" effect.

Figure 2 .

 2 Figure 2.19: The leaf geometry from Figure 2.13e is lifted to 3D, and combined with a set of randomly placed spheres through the hydrophobic contact operator from Figure 2.16.

Figure 2 .

 2 20a. This situation arises when different sketches define constraints that are conflicting in the operator domain. For instance, in Figure 2.20a, the top left and top right fields are mirrored by vertical symmetry. In that case, pairs of symmetric points in the Euclidean space define the same point in the operator space. If two different sketches are defined in each side, as the sharp angle in the top left and the blend in the top right, the result of our

Figure 2 .

 2 Figure 2.20: (a) The sketches result in conflicting samples in the operator domain D, revealing a limitation in expressiveness of gradient-based operators; the user can nonetheless opt for an operator exhibiting either C 0 or C 1 continuity to approximate its input. (b) The smooth nature of the operators is unsuitable to capture dense high-frequency information, which is simply ignored as noise by the optimization. (c) On the other hand, this could be thought of as an advantage, as an artist can draw multiple noisy strokes, and obtain an operator fitting them well.

  recent years, the computer graphics community has invested exceptional efforts in adapting the (non real-time) physical simulation algorithms at the core of cinematic special effects (e.g.: [Ziva Dynamics]) to the realm of (real-time) interactive applications (e.g.: games, AR/VR). Many of these advancements have been possible thanks to a new class of physics solver, pioneered by[START_REF] Müller | Position based dynamics[END_REF], realized on top of Verlet-class integrators[START_REF] Bender | Position-based simulation methods in computer graphics[END_REF]. These position-based solvers are capable of elegantly modeling constrained Newtonian dynamics, including rigid-body, cloth, ropes, rods and fluids in a unified framework. A primary example is the NVIDIA FLEX system[START_REF] Macklin | Unified particle physics for real-time applications[END_REF], capable of modeling complex and varied physical phenomena in real-time by leveraging modern GPU hardware.

Figure 3 . 2 :

 32 Figure 3.2: (top) Skeletal striated muscle as a collections of nested fascicles, fibres, and myofibrils; base image courtesy of [Lee et al., 2010]. (bottom) We abstract the fascicles as a collection of rods. These can be overlapping, and their rest-pose structure is controlled by shape-matching constraints.

Figure 3 . 3 :

 33 Figure 3.3: Bio-mechanically accurate simulation of volumetric anatomical structure is the most effective way to simulate secondary motions (e.g. skin sliding on muscles) and deliver true realistic appearance to dynamic virtual characters; image courtesy of [Ziva Dynamics].

Figure 3 . 4 :

 34 Figure 3.4: (top) The parameterization of a VIPER rod, and its discretization. (left) Its rest configuration, and (right) a deformed configuration.

  attached to every point c(z) ∈ R 3 on the curve. The orthogonal frame D(z) = s(z)R(z) is a combination of a uniform scale s(u), and an orthonormal matrix R(z) = [u(z), v(z), w(z)]; see Figure 3.4. Note that our model generalizes classical elastic rods [Spillmann and Teschner, 2007; Kugelstadt and Schömer, 2016], as in those models the scale is kept constant along the curve. As shown in Figure 3.4, any point in the parametrized volume of the rod can be transformed to the rest configuration by a function p(x, y, z) : R 3 → R 3 :

Figure 3 . 5 :

 35 Figure 3.5: Static behavior -We compare the deformation of a standard cosserat rod to our volumetric invariant version. Our additional degrees of freedom allow us to model the buckling (resp. bulging) caused by the stretching (resp. compression) of the rod.

Figure 3 . 6 :

 36 Figure 3.6: Dynamic behavior -Our generalized rods do not only capture static volumetric deformations, but the scale's velocity allows us to model volume dynamics. In this example, note the rod length is unchanged, but our formulation models a volumetric shockwave travelling through the rod.

  .16) where K b = [ k x b e x k y b e y k z b e z ]. This energy compares the Laplacians of the deformation functions giving a second order measure of the deformation. Volume potential. By denoting the determinant with | • |, and stiffness by k v , we define the volume preservation at midpoints as:

Figure 3

 3 Figure 3.8: (left) An elastic band mesh and its VIPER discretization. (right) The VIPER simulation and their deformation "skinned" to the rest-pose mesh.

Figure 3 . 10 :

 310 Figure 3.10: Muscle contraction -(top) Muscle at rest and its excitation (force shortening edge lengths in an area) with traditional shape-matching constraints (middle) vs our novel scaled shape-matching constraints (bottom).

  Figure 3.11: The mesh-to-VIPER conversion process. Given source/sink constraints, we compute a harmonic function in the volume, and extract a few discrete iso-levels. Within each of these, we execute a restricted CVD to place 5 elements of the same radii on these surfaces. We then execute a combinatorial optimization that connects samples across layers to produce minimal length curves.

Figure 3 .

 3 Figure 3.12: We illustrate several examples of the VIPERs extracted by our automated process. For the pectoral, our VIPER model employs 8 rods, each discretized by 8 elements. In comparison, the surface meshes by ZIVA contain (|∂V | = 1584, |∂F | = 3164) on the boundary and its (volumetric) simulation mesh contains (|V | = 350, |T | = 1089) tetrahedral elements. Where ∂ ≡ "boundary", F ≡ "faces", V ≡ "vertices" and T ≡ "tetrahedra".

Figure 4 . 1 :

 41 Figure 4.1: An example of MIST architecture -A network H η estimates locations and scales of patches encoded in a heatmap h. The locations of its top-k responses are extracted by E K .Patches are then extracted via a sampler S, and then fed to a task-specific network T τ . The results are then fed to a task-specific loss L. In this example, the specific task is to re-synthesize the image as a super-position of (unknown, local) basis functions.

  MIST architectures. Therefore, we introduce a new deep framework which we name Multiple Instance Spatial Transformer or MIST for brevity. From a high level, the MIST framework first decomposes the image into a finite collection of patches, and then processes these patches to perform a given task. As illustrated in Figure4.1 for the image synthesis task, given an image we first compute a heatmap via a deep network whose local maxima correspond to locations of interest. From this heatmap, we gather the parameters of the top-K local maxima, and then extract the corresponding collection of image patches via a re-sampling process. With the collection of patches, we execute the same task-specific network whose output is aggregated to finally evaluate a task-specific loss. We then optimize this task loss to train the entire framework.Training MISTs by lifting. Training a pipeline that includes a non-differentiable selection/gather operation is non-trivial. To solve this problem we to lift the problem to a higher dimensional one by treating the parameters defining the interest points as slack variables, and introduce a hard constraint that they must correspond to the output that the heatmap network gives. This constraint is realized by introducing an auxiliary function that creates a heatmap given set of interest point parameters. We then solve for the relaxed version of this problem, where the hard constraint is turned into a soft one, and the slack variables are also optimized within the training process. Critically, our training strategy allows us to have an optimizable version of 1 non-maximum suppression, and 2 top-K selection, thus creating a network architecture resembling compute strategies that were dominant in pre deep-learning computer vision.

  not need to normalize here, as (4.2) has unit sum along the s dimension. These two operations are abstracted by our top-K extractor {x k } = E K (h) in (4.1).

. 6 )

 6 Algorithm 1 Multi-stage optimization for MISTs Require: K : number of patches to extract, L task : task specific loss, I : input image, η : parameters of the heatmap network, τ : parameters of the task network. 1: function TrainMIST(I, L task )

  order for the network to learn to only respond to the desired locations, we need negative examples just as much as we need positive examples. To circumvent this problem, we propose a multi-stage training optimization.

Figure 4 . 2 :

 42 Figure 4.2: MNIST character synthesis examples for (top) the "easy" single instance setup and(bottom) the hard multi-instance setup. We compare the output of MISTs to grid, channel-wise, Eslami et al.[START_REF] Eslami | Attend, Infer, Repeat: Fast Scene Understating with Generative Models[END_REF] and Zhang et al.[Zhang et al., 2018c].

Figure

  Figure 4.2 (top). Specifically, we perturb the digits with a Gaussian noise centered at each grid

Figure 4 . 4 :

 44 Figure 4.4: Unsupervised inverse rendering of procedural Gabor noise. In the last column we highlight localization mistakes.

Figure 4 .

 4 Figure 4.5: Two qualitative examples for detection and classification on our Multi-MNIST dataset.

  Figure 4.4 illustrates the results of our experiment. Note how the learnt auto-encoder learnt very well to reconstruct

2 A 4 4 e x πr 4 4 e y πr 4 2f

 244 T Adxdy can be simplified to a block diagonal matrix of the [Rq] T × (Rq) = 0. Moreover, as the center of mass of the disc is placed at the origin D mRqdxdy = 0 and D m[Rq] × = 0. We can now integrate the diagonal elements leading to e z R T is the second moment of area of a disc in world-space. ext dxdy is the sum of the external forces which act on the disc, τ ext = s D (Rq) × f ext dxdy is the sum of the external torques and γ ext = D (Rq) • f ext dxdy is a quantity which can

: Volume Invariant Position-based Elastic Rods

1. Sketch-based Implicit Blending B. Angles, M. Tarini, B. Wyvill, L. Barthe, A. Tagliasacchi Transactions on Graphics, 2017 2. VIPERB. Angles, D. Rebain, M. Macklin, B. Wyvill, L. Barthe, J. P. Lewis, J. Von Der Pahlen, S. Izadi, J. Valentin, S. Bouaziz, A. Tagliasacchi In Transactions on Graphics (currently under submission to SIGGRAPH 2019) 3. MIST: Multiple Instance Spatial Transformer Network B. Angles, S. Izadi, A. Tagliasacchi, K. M. Yi In ArXiv (current under submission to CVPR 2019)

  While our MISTs enable us to approach several vision tasks, such as the ones above, with minimal supervision the true challenge lies in the definition of an effective training strategy. Back-propagation through a selection process is possible, alike to what is performed for max-pooling. In Section 4.6, we argue why this is highly detrimental, and propose an effective multi-stage training solution.

	Evaluation and implementation. We qualitative and quantitatively evaluate MISTs on a
	number of tasks in Section 4.7, and provide further implementation details in Section 4.8.

3

the latent coefficients representing each instance. Further details and additional example tasks are described in Section 4.5.

Training MISTs.

Table 4 .

 4 2: Instance level detection and classification performance on the MNIST hard dataset.

		MIST channel-wise
	IOU 50%	84.6%	25.4%
	Classification 95.6%	75.5%
	Both	83.5%	24.8%

From now on, we will refer to "striated skeletal muscle" simply as "muscles".

Acknowledgements xv

Acknowledgements

I would like to express my deep gratitude to my advisors Andrea Tagliasacchi and Loïc Barthe. I am grateful to Andrea Tagliasacchi for the tremendous amount of time he has invested in advising me. I have learnt a lot by working with him and his work ethic has inspired me to give the best of myself during my PhD and forward. I thank Loïc Barthe for believing in my potential by offering me the opportunity to pursue a PhD and setting up this jointly supervised program for me. I am grateful to Brian Wyvill for being my advisor until his retirement. He has been very supportive of my work and has encouraged me to submit to SIGGRAPH Asia, which turned out to be a great piece of advice. I would also like to thank Kwang Moo Yi for his invaluable guidance on my deep learning work and for being always available to patiently answer the numerous questions I came up with. I have been very fortunate to meet great researchers during and before my PhD. I am particularly

Chapter 3

Simulation with Primitives Figure 3.1: (left) We compactly model muscles as a collection of generalized rods, where volume conservation is expressed by a radius function defined on curve's vertices -vis sphere's radii. (middle) In the limit, this discretization represents the smooth rods illustrated here, for which physical constraints due to volume invariance can be expressed analytically. (right) The rods create a subspace on which physics is solved, and its effects later propagated to the muscle mesh via linear blend skinning; please see the animation in our supplemental video. The "Max" anatomical model is courtesy of ZIVA Dynamics.

Abstract

We extend the formulation of position-based rods to include elastic volumetric deformations. We achieve this by introducing an additional degree of freedom per vertex -isotropic scale (and its velocity). Including scale enriches the space of possible deformations, allowing the simulation of volumetric effects, such as a reduction in cross-sectional area when a rod is stretched. We rigorously derive the continuous formulation of its elastic energy potentials, and hence its associated position-based dynamics (PBD) updates to realize this model, enabling the simulation of up to 26000 DOFs at 140 Hz in our GPU implementation. We further show how rods can provide a compact alternative to tetrahedral meshes for the representation of complex muscle deformations, as well as providing a convenient representation for collision detection. This is achieved by modeling a muscle as a bundle of rods, for which we also introduce a technique to automatically convert a muscle surface mesh into a rods-bundle. Finally, we show how rods and/or bundles can be skinned to a surface mesh to drive its deformation, resulting in an alternative to cages for real-time Chapter 4

Learnable Primitives

Abstract

We introduce a deep network that can be trained to tackle image reconstruction and classification problems that involve detection of multiple object instances, without any supervision regarding their whereabouts. The network learns to extract the most significant top-K patches, and feeds these patches to a task-specific network -e.g., auto-encoder or classifier -to solve a domain specific problem. The challenge in training such a network is the non-differentiable top-K selection process. To address this issue, we lift the training optimization problem by treating the result of top-K selection as a slack variable, resulting in a simple, yet effective, multi-stage training.

Our method is able to learn to detect recurrent structures in the training dataset by learning to reconstruct images. It can also learn to localize structures when only knowledge on the occurrence of the object is provided, and in doing so it outperforms the state-of-the-art.

Introduction

The ability to find multiple instances of characteristic entities in a scene is core to many computer vision applications. For example, finding people [START_REF] Sewart | End-to-End People Detection in Crowded Scenes[END_REF]Zhang et al., 2018a], detecting arbitrary number of classes and objects [START_REF] Ren | Faster R-CNN: Towards Real-Time Object Detection with Region Proposal Networks[END_REF][START_REF] He | Mask R-CNN[END_REF][START_REF] Redmon | YOLO 9000: Better, Faster, Stronger[END_REF], and detecting local features [START_REF] Lowe | Distinctive Image Features from Scale-Invariant Keypoints[END_REF][START_REF] Bay | SURF: Speeded Up Robust Features[END_REF] all rely on this ability. In traditional vision pipelines, selecting the top-K responses in a heat-map and using their locations is the typical way to approach the problem [START_REF] Lowe | Distinctive Image Features from Scale-Invariant Keypoints[END_REF][START_REF] Bay | SURF: Speeded Up Robust Features[END_REF][START_REF] Felzenszwalb | Object Detection with Discriminatively Trained Part Based Models[END_REF]. However, due to the non-differentiable nature of this operation, it has not found immediate application in deep learning based solutions.

Circumventing top-K in end-to-end learning. To overcome this challenge, researchers proposed to use grids [START_REF] Redmon | You Only Look Once: Unified, Real-Time Object Detection[END_REF][START_REF] He | Mask R-CNN[END_REF][START_REF] Detone | Superpoint: Self-Supervised Interest Point Detection and Description[END_REF], to simplify the formulation by isolating each instance [START_REF] Yi | LIFT: Learned Invariant Feature Transform[END_REF], or to provide alternative supervision by optimizing over multiple branches [START_REF] Ono | Lf-Net: Learning Local Features from Images[END_REF]. While effective, they do not generalize well outside the application domain for which they were designed. Other formulations, such as the use method to be trained with also examples where various number of total digits exist (images with only 1 digit, 2 digits, etc.). We make a special exception for this method, and populate the training set with all of these cases; 4 we finally compare to the state-of-the-art method by Zhang et al. [Zhang et al., 2018c] that provides a heatmap-based method with channel-wise strategy for unsupervised learning of landmarks.

Results for "MNIST easy". As shown in Figure 4.2 (top) all methods successfully re-synthesize the image, with the exception of [START_REF] Eslami | Attend, Infer, Repeat: Fast Scene Understating with Generative Models[END_REF]. As this method is sequential, with nine digits the sequential implementation simply becomes too difficult to optimize through. Note how this method only learns to describe the scene with a few large regions. Quantitative results are summarized in Table 4.1.

Results for "MNIST hard". As shown in Figure 4.2 (bottom), all methods except ours fail to properly represent the image, where not only it was able to reconstruct the image, but also learnt how to localize the digits. Note that while it might look like the grid method succeeded, its trained auto-encoder simply failed in capturing the concept of individual digits. Conversely, as shown in Figure 4.3, our method is able to learn this, demonstrated by the auto-encoder successfully separating the existing overlaps. For quantitative results, please see Table 4.1.

Finding the basis of a procedural texture. We further demonstrate that our methods can be used to find the basis function of a procedural texture. For this experiment we synthesize textures with procedural Gabor noise [START_REF] Ares Lagae | Procedural noise using sparse gabor convolution[END_REF]. Gabor noise is obtained by convolving oriented Gabor wavelets with a Poisson impulse process. Hence, given exemplars of noise, our Simulation of primitives. In Chapter 3, the use the "tapered capsule" as a useful primitive for simulation was introduced. One might wonder whether there are more primitives that are well suited for simulating physical objects. An interesting example to investigate is the "wedge"

primitive that is the two-dimensional counterpart of the tapered capsule. Its skeleton is a triangle (as opposed to a segment for a capsule).

MIST for unsupervised object detection. Presented in Chapter 4, MISTs are a first step towards the definition of optimizable image-decomposition networks that could be extended to a number of exciting unsupervised learning tasks. Amongst these, we intend to explore the applicability of MISTs to unsupervised detection/localization of objects, facial landmarks, and local feature learning.

Learning three-dimensional primitives. Another interesting venue for research would be to apply MISTs to three-dimensional data. This would have applications in scene compression, scene understanding, object detection, and keypoints detection.

Animation decomposition. So far, in Chapter 4, we have only used fixed images to learn to primitives. But in an animated sequence such as a video or a three-dimensional animation, primitives have a temporal coherency. This could be leveraged to learn more efficiency and solve ambiguous situations. During an animation primitives may also deform, which should also be taken into account.

Appendices

Appendix for Chapter 3

A.1 Elastic Potentials

Strain. As defined in Section 3.5 we can write the strain energy as

where we drop the subscript from K for brevity. This equation can be extended as

where

We can derive the gradient operator for each part of this summation leading to

where Ω = [Ω u Ω v Ω w ] T is the Darboux vector; and analogous expressions can be derived for ∇p.

We can now observe than (A.1) can be broken up into the sum of distinct energies

as the cross terms evaluate to 0. After integrating over the disc we can reformulate (A.9) as is the second moment of the area of a disc scaled by the stiffness.

Volume. As defined in Section 3.5 we can write the volume energy as

Based on the strain derivation, the determinant can be computed as

We can now integrate over the disc leading to

A.2 Prediction Step

As described in (3.9), the inertial potential over the disc is of the form Center update. By solving the linear system (A.26) for ∆c we find the prediction update for the center ∆c = h ċ + h 2 πr 2 m ξ ext , (A.30)

As we integrate on a disc located at a midpoint this update is valid for the center of the disc located at this point. We approximate the update over the end points by using the same update rule.

Scale update. Similarly, the scale update can be computed solving the linear system for ∆s leading to ∆s = h ṡ + 2h 2 πr 4 m γ ext , (A.31)

We also approximate the update over the end points using the same update rule

Frame update. The frame update can be computed by solving the linear system for ∆θ leading to ∆θ = h θ + (A.36) are the inertial predictions for the different degrees of freedom. The variational form of implicit Euler (3.9) can then be written in the form min

T and X are vectors containing all the degree of freedoms and their predictions, K is a block diagonal matrix stacking the stiffness parameters scaled by the length of the piecewise elements, A is a block diagonal matrix stacking the inertia weights scaled by the length of the piecewise elements, and W(X) = [ W1(X) W2(X) ... ] T stacks the potential energy functions.

Variational Solver. To solve this optimization we can linearize the elastic potentials and write an iterative Gauss-Newton optimization min

where k is the iteration number, X k = X k-1 + ∆X, and we initialize X 0 = X. Since Equation A.38 is quadratic in the unknown ∆X, we can minimize it with a single linear solve

However, the conditioning of this linear system is greatly dependent on how stiff are the elastic potentials. Following the optimization trick presented in [START_REF] Ian | On the accurate determination of search directions for simple differentiable penalty functions[END_REF], for elastic potentials with large stiffness a better option is to split the equation above as

Note that when the elastic potentials are infinitively stiff K -1 vanishes λ = [λ 1 , . . . , λ n ] T becomes the vector of Lagrange multipliers. We can now reformulate (A.39) as ∆X = -h 2 A -1 ∇W(X k-1 ) T λ k -(X k-1 -X) (A.41)

≈ -h 2 A -1 ∇W(X k-1 ) T ∆λ, (A.42) by assuming ∇W(X k ) ≈ ∇W(X k-1 ), and where λ k = λ k-1 + ∆λ and we initialize λ 0 = 0.

This can be proven by induction knowing that X 0 = X and λ 0 = 0. By substituting ∆X in the Equation A.40, we can rewrite the system of equations as    ∆X = -h 2 A -1 ∇W(X k-1 ) T ∆λ, (h 2 ∇W(X k-1 ) T 2 A -1 + K -1 )∆λ = W(X k-1 ) -K -1 λ k-1 .

Therefore, ∆X and ∆λ can be found with a single linear solve.

A.4 Closed form pill projection