Gr Égoire Dupont De Dinechin
email: gregoire.dupontdedinechin@mines-paristech.fr

Alexis Paljic
email: alexis.paljic@mines-paristech.fr

Demonstrating COLIBRI VR, an Open-Source Toolkit to Render Real-World Scenes in Virtual Reality

Keywords:

This demonstration showcases an open-source toolkit we developed in the Unity game engine to enable authors to render real-world photographs in virtual reality (VR) with motion parallax and viewdependent highlights. First, we illustrate the toolset's capabilities by using it to display interactive, photorealistic renderings of a museum's mineral collection. Then, we invite audience members to be rendered in VR using our toolkit, thus providing a live, behindthe-scenes look at the process.

INTRODUCTION

Motivations

Many immersive experiences rely on virtual environments and assets created from sets of real-world photographs. Indeed, images and videos are a convenient way of capturing the real world's scenery for VR, e.g. to create immersive virtual visits, high-fidelity training environments, and vivid entertainment experiences. These imagebased VR experiences are made even more compelling when the photography-based environments and assets provide motion parallax -i.e. the sensation that, when moving one's head, closer objects move faster than objects further away -and render specular highlightsi.e. convincingly expose bright reflections on specular surfaces by adapting objects' colors based on the viewpoint from which they are seen (see Fig. 1). Multiple recent works have thus demonstrated how to achieve motion parallax and view-dependent effects for imagebased VR scenes rendered from both 360° [START_REF] Huang | 6-DOF VR videos with a single 360-camera[END_REF][START_REF] Serrano | Motion parallax for 360°RGBD video[END_REF] and perspective [START_REF] Chen | Automatic generation of dynamically relightable virtual objects with consumer-grade depth cameras[END_REF][START_REF] Overbeck | A system for acquiring, processing, and rendering panoramic light field stills for virtual reality[END_REF] image datasets, by estimating underlying geometry and applying different image-based rendering techniques.

However, despite this interest for creating virtual reality content from photography-based assets and environments, there seem to be very few integrated, open-source interfaces dedicated to providing easy access to image-based modeling and rendering tools. Indeed, although 3D reconstruction is now commonly used by VR content creators, e.g. using consumer photogrammetry toolkits, the reconstruction workflow is generally unadapted to rapid prototyping by small research teams, as it is often costly and complex. As for viewdependent rendering, its study and use seems to be mostly confined to an expert segment of the graphics community, and it is only rarely examined in the context of VR user experience research: instead of dynamically updating the reconstructed scene's color from the closest images to render specular highlights, color data is commonly averaged into static texture maps during reconstruction, which is likely to result in lower-fidelity output views at render time. ©2020 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. Because we use a view-dependent rendering method, specular highlights appear on the mineral's facets as the viewpoint changes.

Our approach

To help tackle this issue, we developed the Core Open Lab on Image-Based Rendering Innovation for Virtual Reality (COLIBRI VR), an open-source, integrated resource that provides image-based rendering methods adapted for VR, based on which to discuss, prototype, and lead research related to user experience in photography-based virtual environments. Because the toolset and codebase are openlyavailable, authors wanting to learn more about how to implement image-based modeling and rendering algorithms now have concrete examples to work with, that they can try out on their data and eventually improve. This also means that results obtained using the toolkit can more easily be compared with and replicated by other research teams. Furthermore, because COLIBRI VR is integrated as an interface within the Unity game engine, it can be used as a natural addition to VR content creators' usual development pipeline, to develop interactive experiences that combine synthetic virtual objects and photography-based assets. By linking the toolset to open-source software for 3D reconstruction and mesh simplification such as COLMAP [START_REF] Schönberger | Structure-from-motion revisited[END_REF] and Blender, we also provide an easy way for authors to estimate geometry from sets of images from within Unity.

In terms of architecture, COLIBRI VR essentially consists of C# scripts and a series of ShaderLab (Cg/HLSL) vertex/fragment and compute shaders. The first define an overarching framework for the different processing methods and provide the toolkit with a graphical user interface. With the second, we implement modeling and rendering operations, such as mesh generation from depth maps and view-dependent blending: we thus provide implementations of several core image-based rendering algorithms, such as Unstructured Lumigraph Rendering [START_REF] Buehler | Unstructured lumigraph rendering[END_REF]. Because we target interactive virtual reality, one of the challenges during development was to adapt and optimize these algorithms so as to maintain high framerates, in order not to cause discomfort. For the same reason, we had to regularly check that the implemented shaders and scripts effectively enable immersed users to act on the image-based assets in a variety of ways, e.g. grab them, change them to be larger or smaller, or make them interact with synthetic 3D objects, as we expect such interactions to be useful in many typical image-based VR use cases.

RESEARCH DEMONSTRATION 2.1 An interactive VR mineralogy museum

The first segment of our demonstration consists in inviting audience members to try out an interactive VR visualization of the MINES ParisTech Mineralogy Museum's collection of minerals (see Fig. 2). The content for this segment is pre-processed before the demonstration: we designed a virtual environment in which to display the virtualized minerals, and used COLIBRI VR's interface to transform the input photographs of these minerals into scene representations adapted for real-time view-dependent rendering in VR. When trying out the experience, audience members are thus able to examine the minerals in 3D space by walking around them and picking them up, and can observe the specular reflections captured as light bounced off the minerals' facets during acquisition.

The objective, for the museum's curators, is to display the collection interactively and with high fidelity in an online showcase, so as to give far-away viewers a teaser of the physical visit and encourage them to come in person. The experience may also enable users to more closely examine the displayed pieces: for instance, many of the real-world minerals cannot be touched due to being fragile and costly, and may occasionally be unavailable for viewing; the virtualized minerals therefore provide a complementary learning experience, as they can be taken and observed from up close, no matter the current state or location of their real-world counterpart. This segment is also an opportunity to showcase COLIBRI VR's potential in the context of a specific use case. Indeed, the specular highlights appearing on a mineral's facets are an important visual cue for understanding its visual appearance and material properties: without the ability to render these highlights, many minerals thus could not effectively be rendered as 3D assets in a virtual museum.

Live reconstruction of audience members with viewdependent gaze following

The second segment of our demonstration aims to illustrate a custom view-dependent rendering effect enabled by the toolkit, by which a virtualized model of a person will seem to always be looking towards the user in VR (see Fig. 3). This effect is inspired by the Welcome to Light Fields [START_REF] Overbeck | A system for acquiring, processing, and rendering panoramic light field stills for virtual reality[END_REF] demonstration, and is obtained by taking multiple pictures of a participant from different angles, in which the participant stays still but keeps looking towards the camera as it moves to different viewpoints. Using the graphical user interface, we call COLMAP to reconstruct a 3D model from these photographs, and use COLIBRI VR's implementation of viewdependent rendering methods to render the model. By dynamically selecting the images most relevant for the current viewpoint before blending them to create the output view, we thus select images taken from a viewpoint close to the user's current one, which are therefore images in which the person seems to be looking towards the user. In this way, the rendered person's gaze will therefore seem to follow the immersed users as they move around the virtual space. This segment is an opportunity to answer the audience's questions and present COLIBRI VR's processing and rendering interface in more detail. Because the reconstruction operations are not real-time, we use one computer to render the museum demonstration for the audience, and another to process the photographs and prepare them for rendering. Once the data is ready, we then show the rendered participant in VR inside the head-mounted display, before starting the demonstration anew with a different group of audience members.

Figure 1 :

 1 Figure 1: Close-ups of an amethyst rendered using COLIBRI VR.Because we use a view-dependent rendering method, specular highlights appear on the mineral's facets as the viewpoint changes.

Figure 2 :

 2 Figure 2: Using COLIBRI VR to render minerals as interactive assets in a VR museum. Top: a pyrite in its museum stand, and photographs of the pyrite used as input for the process. Bottom: the rendered 3D mineral and the toolkit's interface in the Unity game engine.

Figure 3 :

 3 Figure 3: Using COLIBRI VR to render custom view-dependent effects, such as having a virtualized character's gaze follow the user.

ACKNOWLEDGMENTS

The authors wish to thank Eloïse Gaillou and Didier Nectoux from the MINES ParisTech Mineralogy Museum for their support on the virtual reality mineralogy museum project.


published or not. The documents may come

An ontology and a memory island to give access to digital literature works

1 An ontology and a memory island to give access to digital literature works Yan Rucar, Jean-Gabriel Ganascia Abstract: we have designed an ontology to index a corpus of digital literature works. We have given this ontology the shape of a memory island, a navigable virtual territory where categories are regions and descriptors places, and where archives of these ephemeral works are made accessible.

A corpus: a technically and conceptually digital literature

To define our concept of digital literature, we will first state what it is not. Not all texts available in an electronic form pertain to the field of digital literature, and distinguishing a native digital literature is no more relevant. As a matter of fact, texts which are written in Word, and published online, can very well be printed without their meaning being altered in the least. Are considered as intrinsically digital the literary forms whose development solely relies on computers. The transfer of these works onto a paper leaf would compromise the integrity of their structures, while weakening their contents.

Digital poetry revolves around words moving on the screen, hypertext fiction relies on hyperlinked pages arranged in a tree structure, the generator builds a text through a computer program, network-based fictions are updated in real time by Web contents.

These works are comprised of different technical components and media, which also convey meaning. Some interfaces are conditional paths to the texts, because going through them requires deciphering a complex system, a phenomenon that Espen J.

Aarseth has named ergodicity [START_REF] Aarseth | Cybertext: perspectives on ergodic literature[END_REF]. Our task consisted in indexing these works woven with various technical functionings and possessing underlying conceptual depths. We have also made up for the scattering of these works which, despite a few specialized journals (BleuOrange, The New river), are mostly displayed on their authors' websites. Indeed, there are but a few digital literature publishing houses: Eastgate for example. Many portals are accessible (NT2 laboratory database, Electronic Literature directory), where hundreds or even thousands of links to digital works are stored. However, despite this praiseworthy profusion, we have chosen to investigate only one trend of digital literature aesthetics, so as to display a coherence that cannot be found anywhere else. Also, one of our research axes is a feature common to all digital literary works: ephemerality. Because of the quick obsolescence of software and disintegration of platforms, most of these works have vanished since their inception in 1952, with the generator Love letters by Christopher Strachey.

The issue of access guides our research:

1. through a unique portal gathering initially scattered works falling within the same aesthetics 2. through an indexation of the works providing a privileged channel to their meaning 3. through the creation of archives preserving a memory of the works despite their predictable disappearance

Indexation: building an ontology

We have chosen to index digital works that challenge reading, as elusive texts are common in electronic literature, a strategy that Philippe Bootz has called the "aesthetics of frustration" [START_REF] Bootz | Vers de nouvelles formes en poésie numérique programmée?[END_REF]. Through fast-paced displays, blurred mobile signifiers, the saturation of words on the screen surface, misleading interfaces or complex paths to the text, the reading of these works includes constant disruptions, sudden suspensions of meaning, often nonproductive questionings, yet eventually leading to a semantic area at some point. The ontology has been built using Protégé [START_REF] Kapoor | A comparative study ontology building tools for semantic web applications[END_REF] and is in the OWL format. The following screenshots provide a thorough illustration of this ontology.

Images: screenshots Ontology 1 to 4

This ontology includes descriptors explaining the type of archive created, the concepts of the indexed works, the functioning of the computer programs, the categories of unreadability or illegibility (the way the text is elusive), literary genres (specifically digital), the mode of interactivity, the existence of other versions in other formats or dissemination methods (the text may be available, albeit formally divergent, in a printed book; the work may have been exhibited as an installation in a gallery), and the media used (textual, sound, visual properties). Descriptors give an account of both the aesthetics and materiality of the works. Therefore the indexation is both a critical analysis and the technical description of the corpus, this latter aspect being fundamental with regards to the archiving process. The specification of the technical facets of the works informs their preservation.

The aesthetical aspect is represented by the categories of concepts (screenshot 1), forms of unreadability (screenshot 2), and literary genres (screenshot 2), which clarify the intentions behind the works. The concept is the core meaning of the work. Although there are numerous descriptors, three conceptual families stand out, from the utter surrender to the whirlpools visible on the screen to the quest for control over the text. The contemplation of the computer's performance, the technological fascination, is the first conceptual feature, through which the spectator watches the generator cut up the text. The descriptor "immersion in language" defines works that surround the reader/spectator/inter-actor with barely visible words that are likely to collide, which may deprive the users from their bird's-eye view on language, so as to turn this source of meaning into a disorientating agent. The opposite trend, the search for control, is well summed up by the concept of ergodic literature coined by Espen J. Aarseth.

Understanding complex systems is a prerequisite to access the text. This reading mode is the literary equivalent of video game. Literary genres may be classified into four groups: visual poetry, animation, text generator and net-art. Visual poetry and animation plunge the reader into language, while interactivity emphasizes this dimension. The text generator uses the machine as a creative writing device, and net-art weaves together web contents taken from different sources. These literary genres exhibit categories of unreadability pertaining to their own textual constructions. Our notion of unreadability is based on the following definition by Bertrand Gervais: "Unreadability is a situation of deficiency: there is a lack which must be overcome" [START_REF] Gervais | Presbytère, hiéroglyphes et dernier mot[END_REF]. Our corpus presents a superficial unreadability of the text which must be addressed through an investigation led by the reader/spectator/inter-actor.

The descriptors of the types of archives (screenshot 1), workings (screenshot 4), modes of interactivity (screenshot 2), media specificities (screenshots 3 and 4), leave the conceptual field to focus on the materiality of screens. Workings may be generative (the computer program builds the text), interactive, or consist of a sequence unfolding (video film, mapping software), or the computer opening to networks (links to online newspapers). The rules of interactivity may be obvious, even stated, or completely hidden, in which case the work transforms before the helpless reader.

Five types of documents are preserved: texts, images, sounds, processes recorded on video, computer programs (source-codes). The textual dimension of the works (screenshot 3) centres around the processes of appropriation and transformation, fragmentation, the pairing of the signified and the graphic signifier. The sound dimension (screenshot 3) may be central or secondary, the descriptors "sound related to the text" and "speech synthesis" apply to digital works regulated by or regulating sound files. The visual dimension (screenshot 4) may be conveyed by the graphic signifier of the text, or be autonomous yet interconnected because images go through the same motions as the words. Iconicity may also be completely independent from language.

Visualization of categories and descriptors: generation of a memory island

The LIP6 laboratory has developed a system generating memory islands from ontologies.

Categories become territories and descriptors places, according to a concepts visualization process reminiscent of the ancient arts of memory. The rhetorician would memorize the stages of his discourse through loci, images and places located in various parts of a palace. The non-linear navigation on the memory island is similar to this journey through loci, as described by Frances A. Yates: "It is essential that the places should form a series […] so that we can start from any locus in the series and move either backwards or forwards from it" [START_REF] Yates | The Art of memory[END_REF]. This memorization aid also proves useful when it comes to conceptualization, because it makes the abstract tangible, as well as what lies in its vicinity. Developed by Bin Yang and Chakib Chekrouni, our memory island shows territories that are areas of meaning transversal to the indexed works, and descriptors (features) belonging to these categories. When clicked on, these descriptors open onto a list of works pertaining to this specific property, with direct links to them. One first sees the whole island map, i.e. a body of literary works with its distinctive features, before going into the details of the corpus. Accessing the meaning of enigmatic works goes together with their exhibition; the cartography is at once an enlightening tool and a portal. This memory island can be accessed on the Internet1 . Below are close-ups of the island describing the aspects of the works. It is possible to export these views and store them on a computer. There are six zoom levels, as well as a tree structure on the left of the screen, which provides another view of the categories and descriptors of the ontology. When the information icon at the bottom of the screen is activated, a text pops up. This text presents the digital literary aesthetics covered by the ontology and the memory island. Clicking on a descriptor opens a window giving explanations on its function, since part of the critical jargon used in the field of digital literature might be misunderstood by the general public. Contextualization is therefore very strong; not only is the ontology itself very visible, but it is also surrounded by numerous clarifications. These explanatory words will soon be complemented with the authors' answers to questionnaires. One third of the authors have agreed to give their comments, so that the voices of creators counterbalance the bias inherent to every indexation. These questionnaires will soon be accessible through a new region of the island. Some of the archives are already accessible through the memory island, and more documents will be added. All archives, whether already online or soon to be added, consist of:

Screenshot Memory

-568 screenshots -10 source-codes -26 sound recordings -32 video recordings -9 texts Beside the fact that they are traces of an ephemeral literature, these documents have an informative aim. Digital works are ever-changing, any return to a previous stage of these reconfigurations is extremely difficult, if not impossible. The archives offer fixed images and sounds, as well as texts which originally escape reading. The teaching of digital literature would greatly benefit from these stable artifacts.

The memory island, a tool to grasp a diversity of meanings

The ontology and its visualization are non-linear practices (because these representations have as many entry points as categories) which aim to give a critical and technical account of a corpus characterized by a multiplicity of perspectives and material aspects.

Therefore, the conceptualization of this literature through cartography meets creative forms defying any synthetic approach. The fixation of these works as archives is obviously their reduction, but it is the only way to preserve a literature written in sand, a memory enabled by the conceptual framework in which these traces are set.

https://ilesdememoire.lip6.fr/beta/Oeuvres/result/resultid.html


Guillaume Cordonnier

Marie-Paule Cani

Bedrich Benes

Jean Braun

Eric Galin

Sculpting Mountains: Interactive Terrain Modeling based on Subsurface Geology

Keywords: Terrains, Mountains, Interactive Design, Geology. !

Most mountain ranges are formed by the compression and folding of colliding tectonic plates. Subduction of one plate causes large-scale asymmetry while their layered composition (or stratigraphy) explains the multi-scale folded strata observed on real terrains. We introduce a novel interactive modeling technique to generate visually plausible, large scale terrains that capture these phenomena. Our method draws on both geological knowledge for consistency and on sculpting systems for user interaction. The user is provided hands-on control on the shape and motion of tectonic plates, represented using a new geologically-inspired model for the Earth crust. The model captures their volume preserving and complex folding behaviors under collision, causing mountains to grow. It generates a volumetric uplift map representing the growth rate of subsurface layers. Erosion and uplift movement are jointly simulated to generate the terrain. The stratigraphy allows us to render folded strata on eroded cliffs. We validated the usability of our sculpting interface through a user study, and compare the visual consistency of the earth crust model with geological simulation results and real terrains.

INTRODUCTION

V IRTUAL terrains are essential to many applications, from real-time flight simulators and games to movies. They are also used in geology by scientists who study interactions between surface processes (such as erosion) and tectonics, by educators, and by scientists targeting public outreach.

While geologists have been gradually advancing the simulation of mountain formation processes, models are either focused on the mechanics of crustal deformation [START_REF] Willett | Mechanical model for the tectonics of doubly vergent compressional orogens[END_REF] or on computing the Earth's surface response to tectonic uplift [START_REF] Braun | A very efficient o (n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution[END_REF]; their purpose is not to generate visually compelling 3D terrains at interactive rates. In contrast, terrain generation methods in Computer Graphics already achieve visually realistic results and recent advances allow for interactivity and user control. However, Computer Graphics methods often fail to capture large-scale landforms commonly observed in nature. One reason for that is that most Computer Graphics methods focus on modeling the surface of terrains, and do not attempt to reproduce the geometry of folds, faults, or the uplift systems that shape mountainous landscapes and which, eventually, produce patterns that are observed at various scales, from eroded cliffs to regularly spaced valleys. Moreover, the lack of stratigraphic information (i.e., information on how the characteristics of rock change with depth) in scanned real terrain databases prevents examplebased terrain synthesis methods from capturing many of • G. Cordonnier We present the first volumetric model for the interactive generation of visually-compelling, large-scale terrains aimed at Computer Graphics applications. We address the long-standing challenge of providing the user with interactive control of large-scale terrains, such as shaping mountain ranges, while maintaining visual consistency with phenomena observed in nature, including the modeling of multi-scale folds. We based our work on a well-known fact in the geological community that mountain ranges are formed by several interdependent phenomena triggered by tectonic forces, which, as we demonstrate here, can be efficiently modeled by using knowledge-based, phenomenological simulations. We coin the term geologicallyinspired simulation for simulations that exploit first-order phenomenological knowledge from geology.

When tectonic plates 1 collide, the Earth's crust is thickened in an incompressible manner, which causes the surface to be uplifted by isostasy: considering that tectonic plates float on a much less viscous and more dense part of the upper mantle like an iceberg floating in the ocean, any variation in crustal thickness must be accompanied by relative vertical motion of the plate's surface. Collision often leads to the formation of faults or dipping shear zones that cut through the entire crust. Within the uplifting region, shortening is accommodated by smaller scale folding. The wavelength of the folds depends on the stratigraphy of the crust, i.e., the thickness of layers previously deposited at the surface of the Earth and their relative competency. The surface uplift leads to fluvial erosion resulting in the formation of dendritic drainage systems and valleys. The rate of river incision is controlled in part by the competency Fig. 1. The user controls the motion of tectonic plates by hand gestures (left). We automatically generate a volumetric terrain model that is folded and eroded (large scale view in the middle). The small scale view of the eroded cliffs (right) shows the folded strata that appears on cliffs.

of rocks, which means that the stratigraphy of the crust strongly influences the shape of river profiles and valley walls and may lead, for example, to the formation of local cliffs or ledges that stand out in the landscape.

Our main contribution is a novel, geologically-inspired deformable model that mimics the Earth's crust under tectonic forces. While we do not claim to reproduce complexity of geological models that incorporate formal mechanical calculations coupled with sophisticated erosion models that require heavy computations [START_REF] Braun | Structural evolution of a threedimensional, finite-width crustal wedge[END_REF], we provide an interactive framework for intuitively controlling the creation of plausible terrains in Computer Graphics. Our framework could be also used for interactive illustration purposes in geology. The geological approximations are justified by the interactivity of our sculpting system and by the targeted result: a geologically-plausible large scale eroded terrain surface that hides most of the higher order tectonic processes.

Our multi-layered (or stratigraphic) model [START_REF] Benes | Layered data representation for visual simulation of terrain erosion[END_REF] of crustal material captures the combined effect of compression and folding. We consider that, at the scale of the crust, rocks behave in a ductile (or viscous) manner and deform by uniform shortening. This is similar to the thin-sheet approximation used in geodynamical models [START_REF] Ellis | Forces driving continental collision: reconciling indentation and mantle subduction tectonics[END_REF], based on the observation that -deep in the Earth's crust-rocks deform by ductile creep at very low strain rate. In the upper part of the crust faulting and folding take place. We use a procedural model to add multi-scale folds triggered by the compression field, based on the specific features of the chosen subsurface stratigraphy. The resulting deformation model is used to compute the uplift of the surface. This uplift is then used in a real-time erosion model, enabling interactive visualization of the resulting sculpting of the surface by erosional processes. The additional stratigraphic details of our model (smaller scale layers) enable us to generate multi-scale folds and produce the final, high quality terrain models used for final rendering.

Our method allows hands-on interactive design of mountain ranges: users can quickly draft and then progressively refine large scale terrains through gestural control, while ensuring that the generated mountain ranges are consistent with first-order geological processes. Figure 1 shows an example of an interactive session where the user sculpts and refines multi-layered earth-crust material as if it were plasticine, using displacement gestures applied using finger interaction on a multi-touch display. The first gestures (Figure 1, left) define the tectonic plates and their relative motion. The resulting collision causes surface uplift and the growth of mountains, where the speed of growth is stored in an uplift map. The combination of the uplift and of an existing fluvial erosion simulation is visualized at interactive rates, enabling the user to refine his/her design. High quality rendering of the resulting terrain at different scales (Figure 1, center and right) shows the effect of a nonuniform input stratigraphy on the shape of valley sides.

We validated our method by comparing our results to simulations from geology and to real terrain data. We show that various phenomena appearing at different scales are correctly captured such as the asymmetric large-scale shape of mountain ranges or the small scale, stratigraphicallycontrolled folds that appear along eroded cliffs. We also conducted a user study with artists and non-artists which confirmed the usability of our sculpting tool.

RELATED WORK

In this section, we discuss algorithms for terrain modeling relevant to our method and then review the use of sculpting metaphors for interactive shape design. For a more comprehensive review of terrain and virtual world modeling techniques, we refer the reader to [START_REF] Natali | Modeling terrains and subsurface geology[END_REF] and [START_REF] Smelik | A survey on procedural modelling for virtual worlds[END_REF].

Procedural generation techniques, such as fractals and noise-based generators [START_REF] Ebert | Texturing and Modeling: A Procedural Approach[END_REF], are efficient multi-resolution models for generating terrains. They hold a prominent role in many applications. Despite a few attempts to force fractal generation to fit user-defined constraints [START_REF] Prusinkiewicz | A fractal model of mountains with rivers[END_REF], [START_REF] Stachniak | An algorithm for automated fractal terrain deformation[END_REF], the common problem of procedural terrain generators is their lack of direct control. In contrast other procedural methods model terrains from a set of 3D feature curves used to specify ridges or river networks [START_REF] Kelley | Terrain simulation using a model of stream erosion[END_REF], [START_REF] Hnaidi | Feature based terrain generation using diffusion equation[END_REF], with, however, little help from the system to improve realism. Génevaux et al. [START_REF] Génevaux | Terrain generation using procedural models based on hydrology[END_REF] introduce a method that builds large terrains by combining river sections and terrains into a single geologically-based model. Another function-based approach [START_REF] Génevaux | Terrain modelling from feature primitives[END_REF] defines terrains as hierarchical construction trees, combining primitives with blending, carving, and warping operators. Guérin et al. [START_REF] Guérin | Sparse representation of terrains for procedural modeling[END_REF] introduce a procedural method to represent terrains as elevation functions defined as a sparse combination of primitives built from linear combinations of landform features stored in a dictionary. These features can be extracted either from real world datasets or procedural primitives, or from any combination of multiple terrain models. A common problem with these methods is their lack of control. If the results do not match the design goals, the user cannot easily refine his/her design and needs to re-launch the generation of the procedural terrain with different parameters.

Crust as viscous material

(Φ, q) → u 0 Crust as layered sheets (Φ, q, L, u 0 ) → u Example-based methods explore the reuse of existing terrains in order to improve realism. Many example-based methods are inspired by texture synthesis; for example the method of [START_REF] Zhou | Terrain synthesis from digital elevation models[END_REF] exploits a user-painted coarse image to control the combination of patches extracted from real terrain data. Terrain samples can also be combined along user-drawn mountain ranges [START_REF] Gain | Parallel, realistic and controllable terrain synthesis[END_REF], or can be edited by resketching the silhouettes of mountains from a first person viewpoint [START_REF] Tasse | First person sketch-based terrain editing[END_REF]. Although these methods improve user control while maintaining small scale features from the original terrain data, they generally do not capture large scale landform features commonly observed in real landscapes and do not claim geological realism.

Erosion simulation approaches can be applied to any terrain model to improve the realism of land-form features. Thermal erosion was introduced by Musgrave et al. [START_REF] Musgrave | The synthesis and rendering of eroded fractal terrains[END_REF], later extended by approaches tackling hydraulic erosion [START_REF] Neidhold | Interactive physically based Fluid and Erosion Simulation[END_REF], [START_REF] Benes | Hydraulic erosion[END_REF], accelerated by using computations on the GPU [START_REF] Mei | Fast hydraulic erosion simulation and visualization on gpu[END_REF], [START_REF] Vanek | Large-scale physicsbased terrain editing using adaptive tiles on the GPU[END_REF], applied to layered terrains [START_REF] Benes | Layered data representation for visual simulation of terrain erosion[END_REF], or controlled in an augmented reality environment [START_REF] Reed | Shaping watersheds exhibit: An interactive, augmented reality sandbox for advancing earth science education[END_REF]. These methods accurately capture small scale terrain features but require a terrain model to be eroded, which can induce large scale inconsistencies. In contrast, the work of [START_REF] Cordonnier | Large Scale Terrain Generation from Tectonic Uplift and Fluvial Erosion[END_REF] proposes a geologically-inspired erosion model that achieves the generation of large scale, plausible terrains by combining fluvial erosion with mountain growth. The user controls the process through the input of an uplift map representing the growth rate of mountains. However, the user is not given any help in designing a plausible uplift map, reflecting the usual features of mountain ranges. In our work, we build on this erosion method and combine it with an intuitive, gesture-based method for designing the uplift map, namely sculpting earth-crust material and allowing mountain ranges to grow through the action of plate tectonics.

Plate tectonics simulations has been performed and studied for a long time in the geological community and many sophisticated thermo-mechanical models were proposed [START_REF] Willett | Mechanical model for the tectonics of doubly vergent compressional orogens[END_REF], [START_REF] Braun | Structural evolution of a threedimensional, finite-width crustal wedge[END_REF]. These models require, however, intensive computations. We are not aware of any existing model enabling an interactive 3D simulation of these phenomena. To the best of our knowledge, interactive 3D modeling of plate tectonics was only used in [START_REF] Michel | Generation of Folded Terrains from Simple Vector Maps[END_REF] to generate a procedural terrain from a 2D vector map of mountain summits. There, the terrain is divided into tectonic plates by using Voronoi-based tessellation and an image processing method is used to generate simple multi-scale folds along plate borders, without taking any layered earth-crust model into account. In contrast, we introduce the first 3D interactive simulation of plate tectonics, thanks to a deformable approach that couples a viscous model with multi-scale folding of multi-layered rocks. Our model provides a more accurate simulation of plate collision and folding processes.

Sculpting has long been one of the main expressive modeling metaphors [START_REF] Galyean | Sculpting: An interactive volumetric modeling technique[END_REF]. In addition to physically-based virtual plasticine [START_REF] Dewaele | Interactive global and Local Deformations for Virtual Clay[END_REF], sculpting systems have proven to be successful for editing as rigid as possible objects [START_REF] Sorkine | As-rigid-as-possible surface modeling[END_REF] or shapes with nested features [START_REF] Stanculescu | Sculpting multi-dimensional nested structures[END_REF], and structured shapes, such as architectural models [START_REF] Milliez | Mutable elastic models for sculpting structured shapes[END_REF]. Peytavie et al. [START_REF] Peytavie | Arches: a framework for modeling complex terrains[END_REF] introduced a hybrid model combining material piles and an implicit representation to interactively sculpt 3D terrain models featuring overhangs and caves. Their approach does not rely on geology, however. Other approaches also use tactile devices [START_REF] Lopes | Expeditious illustration of layer-cake models on and above a tactile surface[END_REF] or sketches [START_REF] Natali | Rapid modelling of interactive geological illustrations with faults and compaction[END_REF] to illustrate sub-surface geology. We present the first application of sculpting to geologically-based earth crust material, enabling users to progressively design mountain ranges through gesturebased interaction on a multi-touch table. As in most sculpting systems, deformation automatically propagates in a plausible way in-between the points the user manipulates; here, through the setting of tectonic, earth-crust motion. An important difference from existing virtual-clay or plasticine material is that the earth crust plates are thin shells. While they deform as viscous material and keep a constant volume under compression, they also fold in compressed regions due to their multi-layered nature.

OVERVIEW

In this section we describe the geological background used as the basis for our approach and the choice of our algorithms and data structures.

Plate tectonics in geology

Continental collision is often the end product of the closure of an oceanic basin driven by subduction of the oceanic plate beneath one of the two continents. When continental plates collide, subduction of the lower part of one plate continues (left plate in Figure 3), while the upper layer (the continental crust) made of lighter, granitic rock resists subduction and experiences compression. This compression is accommodated by the formation of a main crustal-scale fault or shear zone, located above the stable, non-subducting plate, that accumulates the deformation by thrusting [START_REF] Willett | Mechanical model for the tectonics of doubly vergent compressional orogens[END_REF]. In other situations, for example where the crust is hotter, the deformation is more diffuse and the crust deforms as a thin viscous layer [START_REF] Ellis | Forces driving continental collision: reconciling indentation and mantle subduction tectonics[END_REF]. In all cases, this causes uniform uplift of the crustal block and we model this growth by approximating the crust as a viscous material (Section 4).

Simultaneously, this block deforms internally by smaller scale faulting and folding. Due to past tectonic events and/or the deposition of sedimentary rocks at the Earth's surface, the continental crust is commonly composed of a number of horizontal rock layers of different mechanical properties. On top of the basement, made of metamorphic and igneous rock, layers of sedimentary rock are deposited over millions of years, due to the erosion of older mountains or to the deposition of carbonate shells in ancient oceans. When compressed during plate collisions, the difference in mechanical properties between the sedimentary layers makes them fold rather than thicken. We compute the multisale folding from a thin sheet approximation of the crust (Section 5). As a result, mountain ranges are typically asymmetrical, with uplift and folding leading to high topography and steep surface slopes on one side of the mountain where the main structure is located, and more subdued elevations and lower slopes on the other side (i.e., the right of Figure 3). The fact that continental crust is made of a number of different rock layers causes different fold scales to appear (top left of Figure 3). The surface slopes resulting from these complex patterns of surface uplift lead to concomitant surface erosion controlled by the rate of fluvial incision [START_REF] Braun | A very efficient o (n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution[END_REF] which, in turn, results from the shear stress imposed by debris transported by fast flowing rivers on the bedrock (Section 6).

Crust

Geologically-inspired interactive simulation

Our goal is to enable the sculpting of large scale terrains (hundreds of kilometers) at interactive rates, where the user's time corresponds to millions of years in real temporal domain. The corresponding simulated time is around 5 million years in our examples.

User input gestures on a multi-touch input device define the shape of tectonic plates and put them in motion (see Figure 2 and the accompanying video). The user first places a reference hand on the table to define the first tectonic plate, serving as an anchor and intended to be stationary. He/she then moves the second hand in contact with the table to define and set the motion of the second plate (defined as a Voronoi region associated with the second hand) and can repeat this gesture to add new plates. Each time a plate is created, crust material is considered to move rigidly in the convex hull of finger positions, and finger motion is interpolated to define a single, rigid transformation, interpreted as the rate at which the plate moves per time unit (Section 4). Assuming that, at this scale, i.e., when geological temporal scales are remapped into a few seconds, the Earth's crust deforms like a soft, viscous material and inertial forces can be neglected, this motion is unaffected by the collision process and can therefore be assumed to remain constant over time -except when the user applies further edits. The output is a transport field denoted as Φ.

Elevation Layers L i

Sub-layers Subsurface stratigraphy (Figure 4) is defined prior to the simulation using a multi-layered, volumetric model of the Earth's crust: n layers of rock L = {L i }, i = 1, . . . , n are defined with corresponding thicknesses 0.1 ≤ τ i ≤ 5 km and viscosity µ i (a measure of the competence of the rock). They represent variations in rock type and mechanical strength at a variety of scales. Thinner layers are used to represent sedimentary material located above thicker, basal rock layers that represent the layering of the crust, for example separating it into a quartz-rich upper crust and a feldspar-rich lower crust. In our implementation, layer thicknesses are spatially uniform, although extending our method to have them vary locally would be straightforward.

≠ µ,τ ≠ k,θ max
To increase visual realism, each layer can be subdivided into 1 -10 sub-layers of varying resistance against erosion. This represents the fractal nature of natural stratigraphy where large competence contrast is likely across thick layers that are themselves composed of thinner layers of reduced competence contrast.

Geologically-inspired simulation consists of a simulation loop that uses three consecutive actions to compute the terrain deformation at each time step (Figure 2). 1) From the user-defined transport field Φ, we compute the elevation rate (or uplift) u 0 of the terrain caused by the compression between tectonic plates, assuming that the earth crust behaves as an incompressible viscous material in the considered geological time step (Section 4). 2) We then compute the uplift field u i of each layer of rock according to the associated multi-scale folding processes. This is performed using an efficient procedural method based on results from the theory of folding of layered materials (Section 5). 3) We apply fluvial erosion and a landsliding algorithm [START_REF] Cordonnier | Large Scale Terrain Generation from Tectonic Uplift and Fluvial Erosion[END_REF] to shape the terrain during mountain growth due to uplift (Section 6). This results in a time dependent surface elevation z(p, t), in which we keep track of the original rock layer at each node position, for rendering purposes.

Time steps: Recall that we are sculpting animated material: the mountains grow and are eroded while they are sculpted. Therefore, there are two temporal loops: the geological simulation loop based on the current user input (which gives us a stationary compression field and thus stationary uplift and folding, then combined with erosion), and the user interaction loop (Figure 2). The compression field, the folds, and then the uplift need to be recomputed only when a new interaction takes place (e.g., when the motion of the plates is modified). Even when not recomputed, they are used to update the terrain elevation at each frame. The geological time simulation step (Section 4) is not the same as the timestep for any of these loops: it is a theoretical time step used for computing the stationary compression field from the current plates motions. Its value is not important since it expresses a stationary situation.

EARTH CRUST AS A VISCOUS MATERIAL

Shortening of the crust is the first phenomenon to be accounted for when tectonic plates collide. The crust at this scale can be assumed to behave as a viscous material, resulting in thickness changes and uplift. To capture this phenomenon, we use a 2D physical simulation that computes a transport field Φ(p) in the region of influence between plates. The transport field is a map Φ : R 2 → R 2 that describes where each point moves subject to the tectonic forces. It is coupled with a geometric volume preservation method to generate the uplift.

1. Regions A i 2. Influence regions Fig. 5. Gestural control: Three moving plates are created and refined by using successive hand motions and their rigid part are set to the convex hull of the interacting fingers [START_REF] Willett | Mechanical model for the tectonics of doubly vergent compressional orogens[END_REF]. The region of influence of each plate is defined as the part of the plane that is closest to its rigid part (2).

Moving plates creation

Every plate is characterized by its supporting region A i where the earth crust displacement is rigid. The associated transformation T i (a rotation and a translation) represents the rigid displacement of the underlying crust material during a simulation step. While A i does not move over time and T i is constant, the continuous displacement of crust material causes surface uplift and mountain growth in colliding regions.

To define a new plate, the user moves his/her hand in contact with the touch table. We provide visual feedback by displaying an arrow from the initial position c i k of each finger to its final position, i.e., when contact with the table ends. The support region A i is defined as the convex hull of the initial finger positions c i k . Without loss of generality, the vector Φ(c i k ) -c i k , corresponding to each arrow, is interpreted as the motion of the underlying earth crust during a given geological simulation step ∆t.

Since the user gesture may not be fully rigid, the plate's rigid transformation T i is computed from the user defined set of transport directions at control points Φ(c i k ) by using least squares minimization:

argmin k Φ(c i k ) -T i c i k 2 .
The term Φ(c i k ), and the corresponding arrows are then updated from T i so that the arrows always depict the rigid motion of the crust material under the rigid part of the plate during a geological step ∆t. The user can use extra gestures to add new arrows and edit T i accordingly.

When crust material moves away from the regions A i to collision regions, the assumption of rigid motion is no longer true. We call the influence region of each plate the portion of the crust closest to the associated A i . We interpolate the displacements T i to compute the transport field Φ in those regions of influence.

Viscous compression

The transport field Φ describes the local movement of earth crust material during ∆t under the action of tectonic forces, and is computed as the interpolation of the transformations T i given by the user. At the space and time scale of our simulation, the earth crust can be approximated as a thin plate of viscous material that shortens under the action of plate motion and collision (Figure 5 (b)) [START_REF] Ellis | Forces driving continental collision: reconciling indentation and mantle subduction tectonics[END_REF]. At this scale rock deformation is isochoric and the viscous deformation can be assumed to be incompressible such that horizontal shortening yields thickness changes.

The transport field. To simulate the 2D compression, we adapt an existing model of viscous virtual plasticine [START_REF] Dewaele | Interactive global and Local Deformations for Virtual Clay[END_REF]. The amount of viscous crust material q is stored in a 2D regular grid, to be transported from cell to cell through Φ. The grid is defined by a step dx and an orthogonal frame (e x , e y ). Each grid cell is represented by its center x i , and q i represents the volume of crust rock inside the cell. The area of the grid cell is denoted by a cell .

A 2 A 3 π(p,A 1 ) d(p,A 1 ) p A 1
Fig. 6. Transport at p is interpolated between support regions A 1 , A 2 and A 3 , depending on the projection π(p, A i ) of p onto T i and distance d(p, A i ) between p and T i .

In our implementation, computations are performed on a 1024 × 1024 resolution grid. Directly implementing a forward displacement approach yields unnatural patterns similar to Moiré'. Therefore, we solve this problem by computing the inverse of the displacement field Φ -1 for transporting material. For each grid point p inside the transformed rigid area T i (A i ) the inverse of the transport is set to the rigid transformation

Φ -1 (p) = T -1 i p. (1)
Outside the region A i , the field Φ -1 is computed by using an interpolation based on a distance weighting operator. Let d(p, A i ) denote the distance between p and region A i and π(p, A i ) the projection of p onto the region A i (Figure 6). We define:

Φ -1 (p) = p + i T -1 i π(p, A i ) -π(p, A i ) d(p, A i ) e / i 1 d(p, A i ) e (2)
and we set the exponent e = 0.9 to get a smooth falloff of the influence of the distance.

Crust material transportation

We use Φ -1 to compute the amount of crust material q in grid cells after the geological step ∆t, i.e., when plate material is transported from the tip to the end of the arrows (Figure 5). Let q 0 be the initial amount of material in each grid cell. In order to obtain the amount of material q(x i ) in a cell c i at the end of the plate's motion, we use the inverse transformation Φ -1 to compute the original positions of the four corners of c i before motion, which defines a quadrangle of area a i . Then we compute the value of q in the cell c i after the simulation step: q(x i ) = q 0 a i /a cell . The initial amount of material is factored later to compute the uplift, so its value does not have any physical meaning. We simply set q 0 = 1.

Uplift from thickness Changes

Crustal deformation is isochoric, i.e., crust material retains a constant volume over time. The local thickness of the crust is proportional to the volume of crust rock inside the associated cell and to the amount of material q in the 2D grid. Therefore, variations in the amount of material provide a map of thickness change rates: dq/dt = (q -q 0 )/∆t.

Since our simulation method does not capture the main fault formation process and the resulting discontinuity in altitude (Section 3.1), we model them by a procedural approach. We define the fault as an offset of the boundary of the compressed region on the stable plate side that is computed as an iso-level of q. then define the uplift due to compression as u 0 = 0 on the stable side of the fault and u 0 = α dq/dt elsewhere. The coefficient α is a scaling coefficient used to set the uplift to the appropriate geological time scale; in our system we set α so that the uplift is 5.0 10 -4 m /y on average. This approach captures the asymmetric nature of the earth crust compression. The resulting uplift map u 0 is further improved by incorporating folds, as described next.

Note that this procedure neglects two potentially important aspects of crustal shortening. First, we neglect isostasy by imposing that the compression-induced shortening results in surface uplift only. In most mountain belts, shortening also results in a vertical depression of the base of the crust (and of the plate). Second, we do not prescribe a dip to the main fault. In nature all compressive faults (or thrusts) dip at an angle of a few tens of degrees with respect to the surface. The consequence of our approximation is that we do not accurately compute the additional horizontal velocity field that arises from the deformation process. We note, however, that such an approximation is also commonly used in the geological community as the effects of horizontal advection of landforms have only been recently investigated [START_REF] Castelltort | River drainage patterns in the New Zealand Alps primarily controlled by plate tectonic strain[END_REF].

EARTH CRUST AS LAYERED SHEETS

The deformation and shape of mountain ranges are also controlled by folding of layers within the crust under compression. More precisely, the crust is composed of several layers of rocks of varying thickness and strength. Complex multi-scale folding occurs, where small folds appear on top of larger ones. Since simulating this complex folding process would be computationally demanding, we use a procedural model to capture this well understood behavior and then rectify the uplift map u 0 , which only captures the viscous behavior of the crust.

Folding of layered materials

Folding is a mechanical instability that results from the growth of stress imbalance at the interface between layers with different competence. Researchers in both geology [START_REF] Yamato | Dynamic constraints on the crustal-scale rheology of the zagros fold belt, iran[END_REF] and Computer Graphics [START_REF] Rémillard | Embedded thin shells for wrinkle simulation[END_REF] have studied this phenomenon. They rely on the observation that the wavelength of folds λ of a layer k of thickness τ k and viscosity µ embedded in a viscous material are correlated, which can be expressed by the following equation [START_REF] Schmalholz | Finite amplitude folding: Transition from exponential to layer length controlled growth[END_REF]:

λ k ∼ τ k τ k-1 1/2 µ k /µ k-1 1/6 . (3)
Our work relies on this theoretical law as it provides a good approximation of the fold formation process. Similar expressions would be obtained assuming an elastic or viscoelastic rheology for the layers; the important point is that we make the wavelength of the folding instability proportional to the thickness of the layer(s) and their competence contrast (viscosity or elastic properties), both of which can be controlled as input parameters.

The user can define the wavelength of the folds that directly controls the folding effects. The fold thickness is computed automatically.

Procedural fold generation

Although folds of appropriate wavelengths could be set to emerge from layered sheet simulation [START_REF] Rémillard | Embedded thin shells for wrinkle simulation[END_REF] or could be generated by combining elastic thin sheet material [START_REF] Narain | Adaptive anisotropic remeshing for cloth simulation[END_REF] with smaller scale folds [START_REF] Rohmer | Animation wrinkling: augmenting coarse cloth simulations with realistic-looking wrinkles[END_REF], we generate the folds procedurally on top of the viscous simulation (Section 4). This allows us to achieve interactive rates, which is key for seamless content creation with sculpting gestures. Moreover, directly reusing wavelengths from Equation ( 3) is sufficient to ensure visually plausible folds. Therefore, we introduce a procedural modeling method for multi-scale folds.

Our method forces the folds to follow a set of lines orthogonal to the compression direction. We constrain those lines to be separated by an average distance λ that is the theoretical wavelength. The lines are chosen based on the compression direction extracted from the transport field Φ. This choice corresponds to the places of maximum folding (called top skeleton below), between which the folds are interpolated. This provides visually plausible results consistent with the expected fold wavelength and height.

Extracting the compression field. We extract the folding direction and amplitude from the crust compression field C by first converting the inverse transport Φ -1 given by Equations ( 1) and ( 2) into a deformation tensor field F:

F = I -∇(Φ -1 -I) -1 .
The stretch tensor field U is computed by taking the square root of the Cauchy-Green deformation tensor [START_REF] Rohmer | Animation wrinkling: augmenting coarse cloth simulations with realistic-looking wrinkles[END_REF]. In our case, this can be written as U = √ F T F. For every cell, the 2 × 2 matrix U i can be diagonalized:

U i = ν 1 e 1 e T
1 + ν 2 e 2 e T 2 , Vectors e 1 and e 2 are the unit eigenvectors in the plane and ν 1 ≤ ν 2 with real positive eigenvalues. We compute a compression field C for the crust, defined in every cell i by: C i = (1 -ν 1 ) e 1 if ν 1 < 1 and C i = 0 otherwise. The vector e 1 is the direction of compression and 1 -ν 1 represents its strength. Geological events that would cause some local extension of crust are beyond the scope of this paper and are not considered. Iterative fold generation. We model the folds resulting from the layered structure by processing layers one after another (Figure 7) and summing the resulting elevation rates u i to get the final uplift rate:

u = u 0 + n i=1 u i . (4)
This captures folding effects at different scales, since each u i represents folding at the specific wavelength λ i (Equation (3)) that results from the interaction between the pair of layers L i-1 and L i . Folds are generated using sinusoidal interpolation [START_REF] Biot | Theory of folding of stratified viscoelastic media and its implications in tectonics and orogenesis[END_REF] between skeletal curves computed among the edges of a 2D grid, and then used to compute u i , as described below. Fold skeletons. We compute crest and bottom lines for the folds called the top and the bottom skeletons. Similar to previous fold generation methods [START_REF] Rohmer | Animation wrinkling: augmenting coarse cloth simulations with realistic-looking wrinkles[END_REF], these skeletons are set to be orthogonal, in average, to the compression direction. We also make sure that they cover the whole compression region. This is done as follows:

We first generate a planar mesh M independent from the terrain grid, and chosen such that the average edge length is equal to the desired wavelength λ i from Equation (3). M is defined as a Delaunay triangulation of a Poisson disc sampling of the terrain, with λ i as average edge length.

Then, some edges of M need to be selected as part of the top skeleton S t (crest lines). We use the compression field C to compute S t as follows. Let [a, b] be an edge of M with center m. For every edge surrounded by two faces F 0 and F 1 , we define a dual edge as an edge linking the center of an edge of F 0 to the middle of an edge of F 1 , both being distinct from [a, b].

We define the compression energy e along an edge [a, b] to measure its alignment with the direction of compression:

e([a, b]) = (b -a) b -a • C i (m) ||C i (m)|| .
An edge with low compression energy is more likely to be in the direction of the folds. Therefore, an edge E is added to S t if an only if its folding energy is below one of its four dual edges. It can occur that the three edges of a face of M belong to S t which would cause a singularity We detect these cases and remove the edge of highest energy from S t . The third step is to build the skeleton S b representing the bottom parts of the folds. Figure 8 shows the template configurations that are used to generate S b from S t . For every triangle F ∈ M we count the number of edges n F of F that belong to S t and apply the following classification.

Template bottom skeleton configurations

If n F = 2 we add an edge to S b , linking the center of the edge of F that is not in S t to the center of F . If n F = 1 we add an edge to S b that links the center of the two edges of F that are not in S t . Otherwise we add three edges connecting the center of one edge of F to the center of F . Note that this construction ensures that the edges of S b remain approximately parallel to the main folds direction, ie. to the average direction of neighboring S t edges.

Fold heights.

In the last step, we build the fold surface between the bottom and the top skeletons. Although the average wavelength is near the theoretical value from Equation (3), it may vary among the folds and we need to know the precise wavelength at any point in order to compute its height. Let d(p, S t ) and d(p, S b ) denote the distance between point p and each skeleton. We define the local wavelength as:

λ(p) = 2 (d(p, S b ) + d(p, S t )) ,
and the local height of the fold at time ∆t (see Figure 9) as: where σ i is the folding speed. We choose it as a fraction (between 0.01 and 0.5) of the uplift u 0 . Note that we add one to the usual cosine curve used to shape the fold [START_REF] Biot | Theory of folding of stratified viscoelastic media and its implications in tectonics and orogenesis[END_REF] because we model the folding behavior of a rock layer that slides over some thicker layer below it, so f i (x) should be always kept positive.

f i (p) = σ i ∆t, ||C i (p)|| 2 1 + cos 2 π d(p, S t ) λ i (p) , S t S b d t d b f λ/2
Local fold refinement. This step aims at complementing our procedural fold model with a local, direct editing tool that allows for authoring the ridges and valleys in a mountain range. Local edits are useful controls for designers in Computer Graphics applications such as movies and games. The user can change the direction of the folds by manually selecting edges and can also change the height of the mountain along a fold line. The advantage of this local refinement tool is that it offers direct control over the local ridges and valleys while ensuring to some extent the geological consistency of the results, since the computed fold skeletons serve as a geological guide for the user to maintain the coherency of the generated mountain.

Uplift update from folds

The uplift u i resulting from the folding of the current fold pair of crust layers is defined as the folding speed ∆f :

u i = ∆f = f (t 0 + ∆t) -f (t 0 ) ∆t .
This value is used in Equation ( 4) to compute the total uplift field u, and yields terrains that fold at multiple scales.

In practice, using several different folding scales is not necessary for the interactive visual feedback during a sculpting session, since only large folds are noticeable. The other fold layers are considered only during the offline rendering of the final terrain.

TERRAIN SURFACE GENERATION

The elevation of the terrain is generated from the combined uplift and stratigraphy of the crust by applying fluvial erosion and landsliding processes. In our system, we did not include a representation of glacial erosion; to our knowledge there is currently no numerical method that can interactively simulate this process which involves the computationally demanding simulation of ice thickness and flow [START_REF] Braun | A new surface-processes model combining glacial and fluvial erosion[END_REF]. km 4 20 Fig. 10. Schematic rendering of a geological cross-section of a synthesized mountain range.

Interactive terrain generation

The terrain surface generation algorithm is based on [START_REF] Cordonnier | Large Scale Terrain Generation from Tectonic Uplift and Fluvial Erosion[END_REF] and it combines the uplift (Equation 4) with a fluvial erosion law. We use the stream power law commonly used in geomorphology that states that the rate of fluvial incision is proportional to local discharge (the flux of water in the river) and bedrock slope [START_REF] Howard | Equilibrium and time scales in geomorphology: Application to sand-bed alluvial streams[END_REF]. Discharge is assumed to be proportional to drainage area. Our method for solving this law relies on a triangle irregular network (TIN) that represents the topology of the terrain. The drainage area A(p) and slope s(p) are derived for each point p in the triangulation in order to compute a rate of elevation change following the stream power law:

∂ z(p, t) ∂t = u(p) -k A(p) m s(p) n , (5)
where coefficients k, m, and n are assumed to be constants (we used n = 1, m = 0.5, and k ≈ 5.5 10 -7 y -1 [START_REF] Braun | A very efficient o (n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution[END_REF]). The erosion equation captures the river carving effect. Moreover, the uplift u is uniformly rescaled so that its order of magnitude is close to 5.0 10 -4 m y -1 , which is a geologicallyconsistent value (uplift rate in mountain belts rarely exceed a few millimeters per year) that produces visually plausible results.

We extend the original model by coupling the fluvial erosion with a landslide erosion model that mimics the way cliffs are shaped by landslides. This method is derived from the original thermal erosion [START_REF] Musgrave | The synthesis and rendering of eroded fractal terrains[END_REF] used in Computer Graphics. We check the slope s(p) after each step of fluvial erosion. If the slope s(p) is larger than a corresponding talus slope θ, we update the height using:

∂ z(p, t) ∂t = l (s(p) -tan θ),
The constant l = 10 -2 m y -1 has been set experimentally and describes the speed of the landslide effect. The talus slope is preset to the commonly accepted default value of θ = 30 • , but our model enables modifications for each layer of rock, to make it depend on material strength.

Considering that the slope at p only depends on the height of its lowest neighbor z(q) and on the horizontal distance between them r = pq , we use:

s(p) = z(p) -z(q) r .
We evaluate z using an implicit time scheme that is computed in linear time, without the need to invert a matrix: this is achieved by evaluating the heights in the right order, ie. z(q, t + dt) before z(p, t + dt). 11. The nature of the rock at the surface, of altitude z, is deduced from its position wrt. fold layers of altitudes h 0 , h 1 and h 2 .

Rock layers at the surface

Elevation at p h 2 h 1 h 0 z Fig.
The multi-layered nature of the crust controls not only the folding behavior, but also impacts the erosion process (see Figure 10). An important feature of our subsurface geologically-based simulation is that we can associate a specific rock layer with each visible point of the eroded terrain surface, and use this information to locally alter the erosional properties or to adapt the rendering.

The rock layer index associated to a point at the eroded surface is obtained from the stratigraphy of the original elevation h 0 of this point inside the crust. This elevation cannot be simply defined as the difference between the total height of uplifted material and the total height of eroded material at that point, because uplift contribution of folds should not be considered when they are being eroded. Therefore, we compute the height of each fold layer at the current simulation step, and evaluate the original height of the surface from its position between the two nearest layers (Figure 11). The height h k of the k-th layer surface is computed by summing the heights of the layers below, starting with h 0 = u 0 t, and adding the specific uplift due to folding of the current layer:

h k (p) = h k-1 (p) + u k (p) t,
where t is the current time in the simulation, since the beginning of mountain growth.

We first compute which crustal layer every point of the surface of the terrain p belongs to. This is achieved by computing the vertical displacement of each rock layer, from bottom to top, until the total thickness gets higher than the current surface altitude z(p); thus the index k satisfies h k-1 < z < h k . We then compute the elevation of the point between the two layers:

γ = h k (p) -z(p) h k (p) -h k-1 (p)
.

The height of the point before deformation is:

h 0 (p) = (1 -γ) h k (p, t 0 ) + γ h k-1 (p, t 0 ),
where h i (p, t 0 ) denotes the altitude of the crust layer i before deformation at time t 0 . We derive the layer from the original height of the point. To imitate geological strata, we divide the layers into sub-layers of different strengths, and adjust the erosion parameters accordingly: we use a high erosion constant k and talus landslide slope (θ = 30 • ) for weak rock layers, and a low erosion constant k with a high slope (θ = 50 • ) for hard rock.

IMPLEMENTATION, RESULTS, AND DISCUSSION

We developed our algorithm in C++ and with OpenGL and GLSL APIs for real-time rendering. High quality images were directly rendered by streaming our synthesized terrains to Terragen 3 which was also used to generate terrain surface details. All examples were generated on a desktop computer equipped with an Intel Xeon E5 CPU, clocked at 3.7GHz with 16GB of RAM. We used a Wacom CINTIQ 24HD multitouch table for the input.

Architecture

We separated the user-interaction and the tectonic simulation into two CPU threads. The first thread handles the user-interaction, the visual feedback, and gesture analysis. The viscous material simulation, fold computation and tectonic simulation were performed in the second thread. The separation of tasks into independent threads allowed us to provide authoring in real-time with an average performance over 30Hz.

λ [m] Skeleton edges Skeletons [ms] Height [ms]

7, 000 955 0.9

TABLE 1 Statistics for the different folds: wavelength (in m), number of skeleton edges, and computation time for the skeleton and the elevation (in ms).

The viscous crust compression is computed on a grid of resolution 1024 × 1024. The algorithm (Section 4) is straightforward to implement on GPU because each grid cell can be computed independently. We achieve high processing performance (10ms) with no visible grid or aliasing artifacts. The folding simulation is divided into two parts: the construction of fold skeletons (Section 5) was implemented on the CPU, whereas the fold height computation using the inverse distance transform was developed on the GPU. Because the distances between the skeleton edges are small and regular, we developed a GPU algorithm that renders small rectangles around each line, no larger than the theoretical wavelength, and output the distance to the line at each pixel. The distance at the junction of two lines is defined as the minimum of the results. The performance of this step is detailed in Table 1.

The terrain surface is a mesh with 10, 000 vertices for a 100 × 100km 2 terrain. The generation process is fast enough to guarantee interactive visual feedback. One iteration takes 10 ms on average and a 5 million years simulation takes only 5 iterations. Thus, we can perform these iterations interactively before updating the terrain on screen. The parameters used in our setup are reported in Table 2. Only a few post processing steps are needed to achieve high quality rendering. First, the elevation is recomputed off-line by combining the uplifts maps resulting from our method with erosion simulation [START_REF] Cordonnier | Large Scale Terrain Generation from Tectonic Uplift and Fluvial Erosion[END_REF] with a smaller time step and a larger resolution (1, 000, 000 vertices). A simulation step runs in 4s in average and we need ≈ 100 iterations to obtain a final terrain. Finally, we add noise weighted by a function of the slope and the elevation to obtain small surface displacements that are shaded during the final rendering step (Figure 13).

Qualitative and quantitative results

Sculpting sessions. Figures 1,[START_REF] Braun | A very efficient o (n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution[END_REF][START_REF] Hnaidi | Feature based terrain generation using diffusion equation[END_REF] and the accompanying video demonstrate the interactive local and global control over the entire process which is one of the important contributions of our method. The user mainly controls his/her design by using hand-on interaction on a multitouch table. By placing and moving both hands on the table, the continental plates are set in motion. The arrows provide visual feedback about the plate shape and location and their rigid transformations are defined by using finger gestures. The use of the fold refinement tool is illustrated in Figure 12 (3).

Resulting terrains. Figures 1,[START_REF] Gain | Parallel, realistic and controllable terrain synthesis[END_REF][START_REF] Génevaux | Terrain modelling from feature primitives[END_REF] and the accompanying video show high quality renderings. Notice the overall asymmetric range, the alignment of the sub-ranges orthogonal to the compression direction, and at lower scales, the folded strata appearing along eroded cliffs. Figure 20 shows the effects induced by the layered erosion properties on the landforms of the terrain.

In addition to producing still versions of high quality terrains, our method can also be used to produce animations of the formation of mountain ranges, as shown in the accompanying video, in cross section and then as a full rendering.

Validation and discussion

This paper is co-authored by an expert in geology, Jean Braun, who helped us validating the modeling choices and the results. A detailed validation follows:

Our simulation

Geological simulation

Fig. 15. Comparison between a terrain profile generated by using the compression step only in our method (top) and a geological simulation result, also neglecting folds (bottom, courtesy of J. Braun [START_REF] Braun | Dynamical lagrangian remeshing (dlr): a new algorithm for solving large strain deformation problems and its application to fault-propagation folding[END_REF]). Note the similar nearly uniform uplifts in the central parts.

Figure 15 shows a comparison between results from geology [START_REF] Braun | Dynamical lagrangian remeshing (dlr): a new algorithm for solving large strain deformation problems and its application to fault-propagation folding[END_REF] and a mountain profile created by our algorithm. It validates the fact that adequately captures the vertical displacement due to crust compression. Figure 18 shows the folding of one layer in comparison to a geological simulation result [START_REF] Fernandez | Fold interaction and wavelength selection in 3d models of multilayer detachment folding[END_REF]. We also provide statistics computed on the graph of these folds. They include wavelength histograms and an histogram of angle between the fold direction and the orthogonal to the compression direction (Figure 16), and a comparison of the graph properties (Table 3). We choose to measure three values that are essential to the visual aspect of the folded regions, namely the number of end nodes, the number of triple joints, and the average fold length (between the nodes or triple joints).

Our results Geological simulation [START_REF] Fernandez | Fold interaction and wavelength selection in 3d models of multilayer detachment folding[END_REF] End nodes We also provide a side-by-side comparison of a result generated by our method with real data from USGS in Figure 19. We selected two mountain ranges that form along plate boundaries, such as the main fault of New Zealand's Alps (Figure 19, top) and the Cerro Chorolque in Southern Bolivia (Figure 19, bottom), i.e. along the eastern side of the Andes. Both formed by thrusting and folding of layered rocks under compression. They differ mostly by their width, with the New Zealand Alps being much narrower than the Andes. Deformation in the New Zealand Alps is focused on a single crustal-scale fault (the Alpine Fault) and built out of a single lithology, the Otago Schist, whereas in the Andes, the stratigraphy is more complex and controls the spacing between the crustal-scale and folds. The main difference between real data and our results is caused by glacial erosion, which is not captured by our method, and results in broader, flat-bottom (or U-shaped) valleys. The importance of including folding and a stratigraphic control on erodibility is evidenced in our examples. Where such folding and stratigraphic control is not included (as shown in Figure 17.1), all valleys strike perpendicular to the axis of the mountain. Although such an organization of drainage can be observed locally, the morphology of most mountain ranges at the largest scale is strongly controlled by crustal-scale folding. In some mountain ranges, folding is mainly controlled by a single layer (Figure 17. Limitations Our sculpting method has several limitations. First, the fold model we use is only valid at early folding stages: when a real mountain fold grows too high, the relation with wavelength used in this paper is no longer true, and the fold does not have a sinusoid shape anymore [START_REF] Schmalholz | Finite amplitude folding: Transition from exponential to layer length controlled growth[END_REF]. A volumetric simulation would be necessary to compute more complex behaviors such as fold distortion [START_REF] Kaus | 3d finite amplitude folding: implications for stress evolution during crustal and lithospheric deformation[END_REF], as well as secondary faults within the crust. While such a simulation would be interesting and would capture a wider range of geological situations in the generated folds, it would be also difficult to control.

Synthesized terrain 1. Real terrain

Second, we do not allow for a plate setup that evolves over time, though it is the case in several real mountain ranges. This would be more geologically correct, but it would provide less intuitive manipulation.

Finally, our method does not lend itself for precise editing as it only provides the user with a high level of control over the final result. Although this limitation is alleviated by the possibility to edit fold skeletons, the user cannot author certain details, such as the number of dendritic patterns while keeping the other unchanged as shown in Figure 20.

User Study

We conducted a user study to evaluate the usability of our mountain sculpting tool. It was tested on 11 subjects who had low to medium knowledge of geology. Among them, three were professional CG artists. The participants first performed an informal learning trial in order to understand the tool. After that, they were asked to reproduce real mountain ranges shown in two images: the Caucasus and the Alps. The users had to evaluate the interactivity and the effectiveness of the method, and the quality of their results. The measurements were rated on a 1-4 Likert's scale.

The quantitative part of the user study showed that the participants had a positive experience while working with our sculpting tool. The effectiveness of the tectonic gestures averaged 3.3, the local control tool received a mark of 2.8 and the overall intuitiveness was 3.0. The users rated their own achievement in reproducing the shape of real ranges with an average of 2.6.

The qualitative part of the user study included some positive comments from the subjects regarding the simplicity of the first contact, the enabled creativity, and the realism of the results. An important feature of the editing tool was the use of arrows to symbolize tectonic motion. The automatic update of arrows during interaction so that they always show the same rigid transformation, was pointed out by users: some found this counterintuitive, whereas most found it to be of great help. The interactivity and simplicity of our system has potential applications in education, in particular for teaching geology: in the field, the geometry of a fold may be partly obscured by the way erosion has revealed it to the geologist. Our system can show the fold as texture on a cliff wall after erosion as well as its true 3D geometry, making it a useful tool for training geologists in reconstructing the fold 3D geometry from its exposed geometry.

CONCLUSION

We have presented the first approach enabling users to progressively sculpt large-scale mountain ranges by using hands-on interaction with displacement gestures, as if they are sculpting plasticine. We use two coupled sub-models, one for modeling the local thickening of the earth's crust due to compression of colliding plates, and the other for modeling multi-scale folding behavior. Interleaving these two sub-models results in a complex uplift map simulating the growth rate of mountain ranges between tectonic plates. This enables us to generate visually compelling terrains at interactive rates.

Independent of the interaction, the terrain the user manipulates evolves according to first-order geological laws.

Consequently, our validations show that terrains modeled using our system capture important landform features of mountain ranges, such as the asymmetry in landform resulting from the presence of a crustal scale fault and multi-scale folding of crust layers appearing in eroded regions.

In future work, we would like to take into account other types of tectonic activity such as volcanoes, and other types of erosion phenomena such as glacial erosion. In addition to the application of our sculpting system to interactive terrain design, we would like to apply them as an educational tool for teaching geology. Lastly, the user input could be complemented by a time-line allowing to create easily tunable animations of the formation of mountain ranges.

APPENDIX A-GLOSSARY

Tectonic plates: large pieces of the outermost layer of the Earth that moves as coherent plates driven by convective motions in the underlying mantle. Two types of plate exist: continental plates which have a crust made of relatively lighter granitic rocks and oceanic plates which have a crust made of relatively heavier basaltic rocks. The motion of the plates is facilitated by the presence of a low viscosity layer just beneath the plates called the asthenosphere. On geological time scales, the asthenosphere appears as a quasi inviscid fluid.

Crust: the crust is a the upper most layer of the Earth; it is chemically different from the underlying mantle. Its thickness varies between 30-40 km in continental areas and 5-10 km in oceanic areas.

Isostasy: according to the principle of isostasy, the tectonic plate can be regarded as "floating" at the Earth's surface on a low viscosity denser "fluid" called the asthenosphere. Due to its low viscosity, the asthenosphere cannot sustain large stress gradient. Any variations in crustal thickness must be accompanied by vertical motion of the surface such that the weigh of any column of rock above the asthenosphere remains constant on geological time scales. This explains why mountains are regions of thickened crust while basins are characterized by an unusually thin crust.

Continental collision: as tectonic plates move at the Earth's surface they interact with each other. Where their relative motion is convergent, a mountain belt is created such as in the European Alps; where their relative motion is divergent a rift develops, such as in East Africa. A few plate boundaries are characterized by a strike-slip or mostly oblique relative motion, such as along the San Andreas Fault system.

Stratigraphy: term that is commonly used to describe the near-horizontal layering of rocks in the crust that is created by the process of deposition or sedimentation. The deep crust is also compositionally layered with quartz-rich rocks near the surface and feldspar-rich rocks near its base.

Competence: resistance of rocks to deformation or erosion due to its relative mechanical strength.

Erosion: process by which the shape of the Earth's surface is altered by transport of rocks, most commonly from regions of high elevation to regions of low elevation. This erosion is performed by flowing rivers and sliding glaciers, but also by more local processes, such as soil creep, landsliding, chemical dissolution or eolian erosion.

Fault and faulting: a discontinuity in the Earth's crust that accommodates motion between two adjacent crustal blocks. Faulting takes place at a large variety of scales (from a few mm to hundreds of km).

Fold and folding: mechanical instability that arises when a strong layer embedded in a weaker material is shortened leading to a quasi sinusoidal pattern of deformation. Folds are commonly observed in geological materials, mostly due to the existence of a horizontal strong layering in the Earth's crust.

Fig. 2 .

 2 Fig.2. User gestures define tectonic plate movement that is used to simulate earth crust as if it was made of a viscous material. The compression and thickness change first yields a coarse uplift map u 0 and then a thin layered strata that produces a refined uplift map u. Fluvial erosion then produces the final terrain elevation h. The overall simulation is interactive and allows for real-time authoring by hand gestures.

Fig. 3 .

 3 Fig. 3. Colliding tectonic plates: plate experiencing subduction (left) with uplift and folding of crust material while the lithosphere plunges down, main fault (middle), and stable plate (right).

Fig. 4 .

 4 Fig. 4. The initial stratigraphy: rock layers L i have different mechanical properties, whereas sub-layers differ by erosion parameters only.

FoldsFig. 7 .

 7 Fig. 7. The different fold layers are stacked and their uplift contributions are summed to obtain the final uplift.

Fig. 8 .

 8 Fig. 8. Folds computation starts by choosing edges for the top skeleton St (gray) over a mesh M (dotted gray) in the best fold direction. The bottom skeleton S b (red) is deduced from St following the template (left)

Fig. 9 .

 9 Fig. 9. The local wavelength value is calculated from the distance from top to bottom skeletons, and used to set fold height.

Fig. 12 .

 12 Fig. 12. Successive steps of a sculpting session. The user starts by setting the movements of the plates (1 and 2), before adding local details by choosing folding directions and prescribing specific elevations (3) to get the final result (4).

Fig. 13 .

 13 Fig. 13. Real time rendering (left) and the same terrain with details generated by Terragen (right).

Fig. 14 .

 14 Fig. 14. The user sculpts a mountain range. The results capture the whole range coherency, and massifs are aligned with the fold lines. Eventually, the folds leave marks after erosion.

1 k

 1 Layer folding speed 0.01 × u -0.5 × u u Target average uplift 5.0 10 -4 m y -Erosion constant 2.5 10 -7 -8 10 -7 m -1 θmax Talus angle 30 • -50 • l Landslide speed 10 -2 m y -2

Fig. 16 .

 16 Fig.[START_REF] Zhou | Terrain synthesis from digital elevation models[END_REF]. Folding validation: Comparison of the angles between fold normals and compression direction (left) and the histograms of wavelengths (right) between[START_REF] Fernandez | Fold interaction and wavelength selection in 3d models of multilayer detachment folding[END_REF] (blue) and our results (red).

Fig. 17 .

 17 Fig. 17. Various 100 km terrains obtained from the same overall uplift but with different folds layout: only viscous simulation (1), a single layer (2, wavelength λ = 7 km), two fold layers (2, λ 1 = 10 km, λ 2 = 5 km) and (3, λ 1 = 15 km, λ 2 = 3 km).

km 10 50 Fig. 18 .

 5018 Fig. 18. Folding validation: Comparison between a geological simulation result from [45] (left), and a top view of our results (right), both rendered using the same shading.

Fig. 19 .

 19 Fig. 19. Comparisons between real mountains (left) and our results (right). Top: we compare the main fault with New Zealand's Alps. Bottom: we focus on eroded folds in the Andes.

  2), such as in the Zagros Mountains in Iran which formed by folding of a sedimentary sequence above a weak salt layer. Various results can be obtained by changing the relative size of each fold layer: Figure 17.3 includes two layers of related fold wavelength (5 and 10 km), and can be compared to the real data from the Andes (Figure 19, bottom), while Figure 17.4 shows a mountain formed from a larger difference in fold wavelength (3 and 15 km). At the smaller scale,

Fig. 20 .

 20 Fig. 20. The difference of erosion properties among rock layers induces different effects: cliffs (1), individual small summits (2), or circus (3).

Fig. 21 .

 21 Fig. 21. Comparison of a photograph showing real folds (1) and our synthesized terrain (2).

TABLE 2

 2 Parameters of our system and ranges of values used in the examples.

ACKNOWLEDGMENTS

This work was partially funded by the advanced grant #291184 EXPRESSIVE from the European Research Council (ERC-2011-ADG 20110209) and by the NSF grant #1606396 Haptic-Based Learning Experiences as Cognitive Mediators for Conceptual Understanding and Representational Fluency in Engineering Education. We also thank Intel and Nvidia for providing computer hardware.

Guillaume Cordonnier is a former student of the Ecole Normale Superieur de Lyon, and obtained his master degree from Institut Polytechnique de Grenoble. He is pursuing his PhD degree at Universit é Grenoble Alpes on geology based terrain modeling. His research interests cover procedural modeling and simulation of natural phenomena. Eric Galin is a Professor of Computer Science at the University Lumire Lyon 2 and researcher at the LIRIS laboratory. His research interests in computer graphics focus on procedural modeling of virtual worlds, inverse procedural modeling, simulating natural phenomena and implicit surface modeling.