Mariama Akodad

Jessica Labour

Erika Nogue

Delphine Delseny

Jean-Christophe Macia

Richard Gervasoni

Benoit Lattuca

Nicolas Nagot

François Roubille

Guillaume Cayla

Vascular complications and bleeding after balloon aortic valvuloplasty performed with or without heparin: HEPAVALVE randomized study

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


Clara Sandrini
email: clara.sandrini@paris-valdeseine.archi.fr

Optimistic Suburbia II -Middle-Class Large Housing Complexes

Session 3: Publicness in middle-class large housing complexes as a new way to examine the premises of cultural encounters and social integration Chairs: Dr. Marie Glaser and Prof. Ellen Braae

Social transformations of inhabited spaces in Sofia, Bulgaria

Bulgarian collective housing complexes, built after the Second World War and during the socialist period, present a singular aspect.

The research group on Bulgarian architecture from the Architecture School (ENSA) of Toulouse conducted surveys in Sofia, Varna and Plovdiv between 2011 and 2018. It sought to understand the planning and construction process of these housing complexes, based on archival sources, interviews with political actors and design architects. It also developed field surveys in order to observe the inhabitant's practices and the transformations of the interior and exterior spaces of the dwellings.

The studies carried out on the city of Sofia were the foundation of the entire scientific approach undertaken on the Bulgarian territory. They were conducted with colleagues from the University of Architecture Civil engineering and Geodesy (UACG), and in partnership with the French Institute of Bulgaria, where they were regularly presented.

A singular socialist planning

Located at the base of the Vitosha Mountain, Sofia was the object of successive urban planning after the Second World War, which considerably enlarged the urban territory by organizing the construction of gilo or micro-districts. These new districts, whose design is similar to the notion of neighborhood units, were intended to be autonomous: they have their own administrative centers, facilities and outdoor spaces. Planned and built at different times, they do not all have the same urban and architectural morphology.

The housing units built in these gilo and micro-districts also respond to an adapted vision of the Soviet economy: the "socialist" system put in place after the Second World War has indeed reserved the right for the inhabitants to own their housing, even if the land had been collectivized. The original property deeds even specify that the inhabitants own a percentage of this collective land, without this being spatialized. This percentage differs, however, between the gilo or the micro-districts, depending on the total surface of the development.

These dwellings were designed by Bulgarian architects who, very often, had the opportunity to do part of their studies abroad. USSR was of course the main destination, but from the 1970s, they also came to western Europe. Thus they mixed cultural models and housing typologies. Soon a diversity of collective dwellings stated to emerge, whereas it was unheard in the former Eastern European countries with planned economy originated from Moscow.

In order to reduce construction costs, these dwellings were delivered unfinished. These state buildings provided the roofing and the underground networks. Each flat was also provided with water, gas and electricity supplies, but the interior finishing and the furnishing of the rooms (sanitary facilities, bathroom, and kitchen) was left to the owners, depending on the materials and furniture available.

Mass housing adapted by the residents

Therefore, the inhabitants participated in the interior design of their dwellings as soon as the buildings were delivered. The shortage of materials spread out these fittings over time and led the inhabitants to use recycled materials from the renovation of buildings in the city center. Each flat was in fact individualized as soon as they moved in, even if the interior finishing work took many years.

The outdoor areas, which are collectively owned by the residents, were also delivered unfinished, apart from the areas dedicated to public facilities, street furniture and children's playgrounds, which were delivered and laid out according to models predefined by the State. During the socialist period, the inhabitants were obliged to arrange and maintain these outdoor spaces one day per month. Beyond this obligation, they developed common spaces, again using recycled materials and adapting the facilities to their needs.

Gradually, the inhabitants were also confronted with the maintenance of the common parts of the buildings. Here again, they organized themselves for routine repairs, cleaning. For major repairs, they called on the know-how of each person, by pooling skills, at the level of construction unit groups and buildings. Thus, the status of owners, combined with the heritage of Ottoman law and the socialist rules of land management, committed them to a process of collective management of the common spaces, both inside and outside.

The dwellings have been subject to visible external renovation since the beginning of their construction. The inhabitants have changed the windows, buying or using reused materials. Some have turned their balconies into an extra room. Others have repainted their flaking facades. And, in the context of energy renovations, many dwellings have been covered with insulation from the outside. This process doesn't apply to the whole building or stairwell. The installation of the insulation, as well as its thickness, depends on the economic resources of the inhabitants.

Therefore, each dwelling is individualized within the housing estate, with more or less extensive transformations to the façade. The micro-districts of the city of Sofia have different aspects, linked to their position in the city, their construction dates and the individual and collective skills of the inhabitants. The proximity to the city center or the existence of a public transport network (tramway, trolley bus or metro) facilitates the movement of reused materials from the city center to the periphery. The microdistricts built later, in the 1980s, are only just beginning to undergo individual transformations: better designed, according to the architects, they do not require energy rehabilitation. But some of Sofia's micro-districts offer different scales of common spaces, both in the rehabilitation of the buildings and in the design of the external spaces.

A hierarchy of Common spaces

Three micro-districts are representative of the way people share collective spaces: Lozenets, to the south of the city center, Drujba, to the east of the city, and Lyulin, located to the west.

The micro-district of Lozenets is part of the continuity of the historical urban fabric: it shapes the urban fabric of a suburb where several periods of construction are represented. This urban fabric is made up of individual houses, apartment buildings and facilities: it appears as a green district with many parks and gardens, a forest crossed by the tramway and the city zoo. In the 17th century, Lozenets was a small village on the outskirts of Sofia, surrounded by forests. It had many Catholic and Orthodox temples, the remains of which still exist at the edge of the district. Deforestation began at the beginning of the 19th century and, gradually, individual houses were replaced by multi-storey brick and concrete buildings. After the Second World War, the infrastructure was built and allowed better communication with the neighboring districts and the old center. And the heart of the district was restructured by one of the first micro-district of the city. Today, numerous developments are appearing in the areas between two construction unit blocks and the ground floor of the buildings are being occupied. Several parts of the exterior areas have been privatized by the inhabitants. This privatization is marked by painted wooden fences which enclose a space that is generally furnished with self-built garden furniture, well-tended plantations and a floor covering borrowed from town center public spaces. Sometimes the whole collective space is open to the view of the passerby, with hierarchical arrangements, as found out around the "VIP club". The outdoor space is not really delimited, but it is composed of vegetable gardens, storage spaces and rest areas. A do-it-yourself path leads to the central shelter, made by the inhabitants with recycled materials. It hosts several social groups depending on the time of day: women in the morning, men in the early afternoon, children in the late afternoon and young people in the evening.

The Drujba 1 micro-district is part of the gilo-district of the same name, and was built in the east of the city center. It consists of 170 housing complexes built between 1965 and 1984 and it is surrounded by major roads, including the road from Sofia to Plovdiv. It shows a diversity of buildings, whose design and construction were spread over a long period of time: several typologies of buildings and housing shape the territory, from the artificial lake to the large boulevards. As in Lozenets, the outdoor spaces are developed, maintained and managed by the inhabitants. The distribution of the different types of buildings according to their functions, and the very location of the residential buildings in relation to each other, create a more or less intimate spaces. Inhabitants' occupations are systematically located close to the buildings: they are concentrated in the grassy areas in-between construction units, which stand as an extension of the dwellings.

To be more precise, these occupations of space are organized by stairwell. The inhabitants informally elect a "stairwell manager" who manages the budget, maintenance and repairs. A local tax is collected by the town hall of the gilo-district and is officially redistributed to each building and then informally to each stairwell. The inhabitants who have more financial means add a share to these taxes which can represent a large amount of the fund raised. This extra money is reflected by the quality of the renovations and the exterior fittings.

The diversity of architectural adaptations is also visible on the facades, where patchworks of coverings, with disparate paints and thicknesses are randomly distributed according to the position of the flats. The transformation of the original balconies is also a common phenomenon: in order to increase the interior surface of the dwelling, many inhabitants choose to fill in their balcony, transforming it into a veranda, a loggia, or an additional room. This proliferation of individual and collective interventions thus organizes a heterogeneous urban landscape. The hierarchy of the occupation of those shared spaces, is structured from the family unit to the group of houses.

The Lyulin gilo-district, in the west of the city center, is the largest in the Bulgarian capital. In the late 1960s, the Sofia Design Bureau (Sofprojekt) received a state order to build a socialist housing complex for 120,000 inhabitants on the site of the small village of Lyulin. The master plan was completed in 1969 and provided for the organization of the gilo-district into eleven micro-districts distributed along a cardo and a decumanus and a "green corridor", linking the "modern center" to the Zapaden park located to the south. This gilo-district was partly built in the seventies and eighties. Constructions kept built after the fall of the socialist regime, but in a haphazard manner and were often left to private investors.

Fig.3 -Transformations of a residential building in Lyulin, photograph by Clara Sandrini

The studies carried out in the residential part of the micro-district number 6 of Lyulin have highlighted an attitude of "laissez-faire" of the public authorities. The project was not carried out as a whole: the buildings that were supposed to face the green space connecting Zapaden Park to the modern center were never built, creating thresholds between residential and green space filled by wastelands and "modern ruins". The green corridor is only a dream plan on paper, as it is now occupied by a landfill, a private car park, a shopping center and slums next to empty luxury buildings.

Within the residential spaces, however, the transformations of public space, common space and private space are taken over by the inhabitants, individually or collectively. The consensus between residents allows them to pool resources to renovate spaces. Tacit rules are created at different levels, as in Drujba : those of the landing, the stairwell, the building or the space in-between the construction units. Several inhabitants can thus establish a common project to install new insulation on the façade, others refurbish a cellar on the ground floor to set up a café while others repair the benches at the entrance to the building.

This care given to the spaces allows the inhabitants to overcome the abandonment of the places by the public authorities since the fall of the socialist regime. It testifies the involvement of the inhabitants in the development of the space and a form of participation, posteriori to the construction of the apartment buildings.

A division of space

The uniqueness of the city of Sofia lies in its cultural heritage, characterized by a heterotopia linked to communist planning and the legacy of Ottoman rules with a particular system of social organizations.

The appropriation and the ephemeral, temporary or durable occupations of collective spaces create the conditions for sharing space and characterize it as the product of an interaction between the collectives of actors who conceived the city, those who built it and those who now occupy it.

The implementation of models throughout the territory according to communist planning, as well as the system of distribution of housing and land, have structured an urban space that seems uniform everywhere. However, each of the districts has a typological diversity which suggests that the Soviet model transposed to Bulgaria has been appropriated by the political authorities and by Bulgarian architects.

In addition, this uniformity is matched by the appropriation of different spaces, from the most collective to the most private, which implies coordination in the implementation of occupations. The social practices indicate an appropriation of individual and collective space linked to a system of ownership of the habitat which includes part of the exterior spaces. The exteriors are thus divided up in space and time of day to allow all the inhabitant groups to use them. They are punctuated by ephemeral, temporary or durable facilities, the construction of which is part of a process of reuse on the scale of the city, and the layout of which shows several levels of privatization.

The result is a social logic: the sofiote space, private, common and public, is thus the representation of a superposition of the coordinated actions of the different collectives of actors opposing the space conceived, planned and built, to the space lived, diverted and occupied. The collective housing complexes thus present a participatory aesthetic which articulates state constructions and private interventions. They show a common space, product of a collective hijacking of the designed space which characterizes an a posteriori participation of the construction, and which shows the different levels of coordination of the users and inhabitants.

Short Bibliography

Fig. 1 -

 1 Fig.1 -Photographic survey of the inhabitants' interventions in Lozenets by Emilie Calvet, Master 1 thesis, ENSA Toulouse

Fig. 2 -

 2 Fig.2 -Building typologies in Drujba 1, by Florie Medus, Master 1 thesis, ENSA Toulouse

  Alexandrova-Guest, Milena. 2006. Habiter Sofia au tournant du XXIè siècle : essai géographique sur l'aménagement des espaces urbains post-socialistes. Phd Thesis. Lyon : Ecole normale supérieure. Jannière, Hélène. 2008. "Planifier le quotidien. Voisinage et unité de voisinage dans la conception des quartiers d'habitation en France(1945-1965)". Strates 14 : 21-38. Momchil, Hristov. 2009. "Biopolitiques de l'habitat socialiste. Le privé comme objet gouvernemental dans les politiques de l'Etat bulgare desannées 1950 et 1960". Histoire@Politique, vol. 7, no. 1 : 4-4.


A Bedel

Y Coudert-Osmont

J Martínez

R I Nishat

S Whitesides

S Lefebvre

Closed space-filling curves with controlled orientation for 3D printing

Keywords: CCS Concepts, Computing methodologies → Shape modeling;, Applied computing → Computer-aided design;

Figure 1: Closed space-filling curves optimized under (a) anisotropy, (b) boundary alignment (c) isotropy inside the E and boundary alignment inside the G. Isotropy and anisotropy are optimized in terms of distribution of angles, while alignment is specified in a vector field (here parallel to the boundary). For the three curves, the number of transitions between the two letters of the logo is minimized. (d) Multi-material 3D printed cube, with isotropic layers optimized by our approach.

Introduction

Many additive manufacturing processes solidify material within each layer along trajectories [ALL * 18]. The way these trajectories cover the layer has direct implications on the fabrication process itself and the object's final properties. In particular, the trajectories impact the quality of extrusion processes such as fused filament fabrication, silicon, ceramic, wire arc, and concrete extrusion. Thus, it is essential to minimize the occurrence of travel moves: motions where the material flow is interrupted while the extrusion device moves from one point to another. Such moves take time and introduce defects where the flow stops and resumes: the print head traveling above the part often leaves traces and spurious deposits.

Therefore, a trend in Additive Manufacturing (AM) focuses on determining trajectories that fill a layer in a single continuous motion [ZGH * 16, DSL18, PBS18, JM20]. In such cases, the trajectory is a cyclic space-filling curve. As layers are now printed in cycles, the start and end positions of the print head can be chosen advantageously (e.g., to minimize travel times), while the extrusion flow is only interrupted where necessary.

Another important criterion is the orientation of the trajectories. These directly impact the surface finish, as the trajectories remain visible and induce anisotropy in surface roughness. Several recent works also note that aligning the trajectories with the principal stress directions improves the mechanical robustsubmitted to COMPUTER GRAPHICS Forum (4/2022). ness [SIM16,TM17,LY17,KWW19a,XLM20,FZZ * 20]. A distinct problem is to attempt to obtain isotropic behaviors, while the deposition tends to follow specific directions (e.g., zigzag and concentric patterns). This applies both to mechanical strength [LXM * 19] and thermal differences during deposition [START_REF] Cheng | Enumeration of additive manufacturing toolpaths using Hamiltonian paths[END_REF].

Finally, another context where precise control over the trajectories is required is when switching materials during a print [PJYP14, XJ18, SMB * 19, SSMVL19]. In such cases, the trajectories have to cover the layer while following a target field, ensuring correct mixtures are achieved despite relatively slow gradients of change, minimizing the number of switches. We refer to this property as trajectory zoning.

In this paper, we consider three different objectives: continuity, orientation, and zoning. The requirements of 1) a smooth, closed, and non-intersecting trajectory, 2) covering the entire surface with everywhere the same spacing while 3) being freely orientable and minimizing crossing between zones forms a multi-objective optimization problem.

Interestingly, a similar challenge has been considered before in the context of Computer Graphics, primarily for interactive artistic creation of maze-like patterns [PS06, WT13]. However, a number of limitations make these techniques less suited for optimizing our objectives automatically across large numbers of slices. We provide further discussion and comparisons in Section 5.2.

We depart from existing methods and approach this challenge by decomposing the optimization into two steps, manipulating different aspects of the problem. The first step is a combinatorial optimizer, exploring Hamiltonian cycles in a specially constructed graph covering the surface. The optimizer starts from an initial cycle and modifies it through operators local to the graph, improving the objective functions. The operators locally update the cycle while preserving its Hamiltonian property. Once the global connectivity is determined, our second geometric optimizer further refines the closed curve's geometry. The result is a smooth, evenly distributed trajectory globally optimized for all objectives.

Our trajectories' orientation can be controlled both in terms of global distribution of orientations (isotropy/anisotropy) and local alignment specified as a vector field. This control may be different in distinct zones. The trajectory can also be optimized to remain in a designated area (zoning) for as long as possible. This objective ensures that multi-material mixtures remain stable while in a specific area.

We consider that a cyclic curve is either a simple closed curve on a surface (the plane) or a cycle in a graph embedded in the surface; the context will make the meaning clear. As the cycles in the graphs become space-filling, the two meanings approach each other.

Related work

Space-filling curves and Hamiltonian cycles in graphs find many uses in Computer Graphics. In particular, they can be used to define coherent linear traversals of graphical data, such as organizing triangle and quad meshes into strips [AHMS96, Tau03, GE04, GLLR11], or improving image compression rates [START_REF] Dafner R | Context-based space filling curves[END_REF].

Space-filling curves have unique aesthetic qualities, which led to their use for artistic creations. Bosch et al. [BH04, Bos10] and Kaplan et al. [KB * 05] explore continuous line drawings using TSP solvers to find a Hamiltonian cycle. Wong and Takahashi [START_REF] Wong | Abstracting images into continuous-line artistic styles[END_REF] stylize images, growing an initial cycle in oriented grids obtained from quad meshing. A similar growth process is performed in [START_REF] Dafner R | Context-based space filling curves[END_REF]. The approach of Pedersen and Singh [PS06] grows space-filling curves from initial simple closed curves (e.g., circles), driving the growth through a variety of user-driven guidance fields. This generates beautiful and intricate patterns. The growth is performed with a simulation scheme akin to a particle system, with repulsion and attraction forces. Sharp and Krane [START_REF] Sharp | Variational surface cutting[END_REF] generate space-filling curves as a side effect of their technique for surface cutting through a similar growth strategy driven by a level-set optimization framework. Growth approaches give control over the final trajectories by locally influencing the process, as demonstrated in [PS06, [START_REF] Wong | Abstracting images into continuous-line artistic styles[END_REF]. Thus, they are a possible option to approach the challenges we seek to address. We discuss the difficulties this presents in Section 5.2.

Graph theory approaches have also been proposed to find covering cycles. These combinatorial methods are very efficient since -on specific graphs -they are supported by a linear time algorithm [Tau03]. This is demonstrated in Akleman et al. [AXG * 13] to generate Hamiltonian cycles along surfaces for various artistic effects. We use a similar approach to initialize our optimizer.

Space-filling curves have properties that make them attractive as deposition trajectories for Additive Manufacturing [JM20]. Chen et al. [START_REF] Cheng | Enumeration of additive manufacturing toolpaths using Hamiltonian paths[END_REF] enumerate all cycles within a regular square grid, evaluating each cycle for the resulting thermal properties during deposition. Lin et al. generate contour parallel fill trajectories that form singly connected spirals within a contour, leading to a continuous fill with smooth trajectories. Soler et al. [START_REF] Soler | A generalized approach to non-layered fused filament fabrication[END_REF] create 3D trajectories as Hamiltonian cycles in grids for wireframe 3D printing. In the context of two-photon lithography, Dahaeck et al. [DSL18] divide a surface to be manufactured in a process akin to a Centroidal Voronoi Tessellation and compute a Hamiltonian path within to obtain a continuous trajectory. Other techniques use space-filling curves to define sparse infills with density gradients (e.g., [START_REF] Kuipers | CrossFill: foam structures with graded density for continuous material extrusion[END_REF]); however, our focus in this work is on dense, solid infills. These techniques do not provide any direct control over the orientation of the generated trajectories.

Method

Our algorithm takes as input a polygon P ⊂ R 2 , a distance 𝑑 that (indirectly) controls the spacing around the curve, an orientation field 𝑈, a set of isotropic zones I 𝑖 ⊂ P, a set of anisotropic zones submitted to COMPUTER GRAPHICS Forum (4/2022). General overview of our approach. We start by covering the input polygon P with a graph and finding an initial cycle. Then, we proceed to optimize the topology of the graph. Finally, we optimize its geometry. Each frame visually represents a step of the method, and is related to a different section.

A 𝑖 ⊂ P, and a partition of P into disjoint material zones Z 𝑖 . P and all the zones are represented as sets of oriented boundary contours and may have holes. P specifies the area to be covered. Without loss of generality, we assume that P is a connected set, as otherwise, we process each connected component independently. The orientation field 𝑈 is piecewise constant. Its underlying data structure is a regular grid 𝑉 covering P, with center vectors having a norm of one -where an alignment is specified -or zero -where no specific alignment is requested. Each vector defines a constant orientation across the square grid cell area. In each of these zones, we seek to obtain respectively a uniform and a non-uniform orientation distribution. Each material zone Z 𝑖 is associated with a material that will be used during the printing process. We assume that two adjacent material zones have different materials (otherwise, they can be merged).

The algorithm outputs a smooth, non-intersecting closed trajectory T covering P. T is composed of 𝑁 vertices (𝑣 𝑣 𝑣 𝑖 ) 1 𝑖 𝑁 and 𝑁 straight line segments between consecutive points. We assume that the trajectory is oriented such that each vertex of index 𝑖 has a successor 𝑖 + and a predecessor 𝑖 -. A per-segment area 𝐴(𝑣 𝑣 𝑣 𝑖 , 𝑣 𝑣 𝑣 𝑖 + ) is defined for each segment; it represents the area around the segment. Multiplied by layer height, it gives the volume to extrude. Dividing 𝐴(𝑣 𝑣 𝑣 𝑖 , 𝑣 𝑣 𝑣 𝑖 + ) area by the segment length gives a notion of the average spacing around.

We optimize the trajectory T to achieve a space-filling property with an interspace close to a nominal value. Intuitively, the trajectory meanders throughout P while covering it uniformly. Formally, this is achieved by minimizing the distance from any point of P to T (coverage) while constraining T to have a specific length. The length indirectly controls the nominal spacing around the curve, as given more length, the curve will get closer to points in P.

We further incorporate the alignment, orientation, and zoning objectives. The alignment objective is specified as non-zero vectors in the field 𝑈 and encourages T to follow these directions where specified. The orientation objective controls the edge orientations' distribution, minimizing or maximizing the distance to a uniform distribution. This optimizes for either isotropy or anisotropy. Zoning minimizes the number of crossings between T and the boundary between adjacent polygons Z 𝑖 ,Z 𝑗 . All these objectives can be freely combined; the orientation objective is independently controlled in every zone. We detail all objectives and their computations in Section 3.1 and the Appendix. An open-source implementation of our approach can be found at https://github.com/ mfx-inria/controllable-space-filling-curve.

Our approach operates in three main steps (see Figure 2):

• (Section 3.2) We build a graph embedded in P, so that the nodes are evenly spaced, have exactly degree three, and so that the edge orientations are biased towards any local alignment specified in the input. • (Section 3.3) After building an initial Hamiltonian cycle in the graph, we perform a stochastic combinatorial optimization of the alignment, orientation, and zoning objectives. • (Section 3.4) We further optimize its geometry for even coverage, smoothness, orientation, and zoning, starting from the resulting space-filling curve while targeting a specific curve length.

submitted to COMPUTER GRAPHICS Forum (4/2022).

Objectives

We formulate our objective functions along the segments of a Hamiltonian cycle T (H-cycle for brevity) covering P.

The objectives we define here are used by both the combinatorial optimizer and the continuous geometry optimizer. They are combined in a weighted sum, always using the same weights.

In the following, E denotes the objective terms, and 𝛾 denotes the weights of the objective terms. We start by giving an intuitive definition of each objective.

• The space-filling objective produces a space-filling curve within the surface of P. It comprises the E CVT and E Len objectives. • E CVT is the Centroidal Voronoi Tessellation objective induced by the segments of T , restricted to the polygon P. It encourages T to be as close as possible to every point of P, maximizing the surface coverage. • E Len controls the length 𝐿 of T . This indirectly controls the spacing between the folds of T , as given more length, the E CVT will redistribute the curve to be close to all points of P. • E Lap encourages T to be smooth; it is a Laplacian term applied to the vertices of T . • E Ali aligns T with the input vector field 𝑈 that specifies alignment vectors (vectors with the unit norm, while zero vectors indicate that no specific alignment is requested). • E Ori controls the distribution of edge orientations, either maximizing (anisotropy) or minimizing (isotropy) the distance to a uniform distribution. • E Zon discourages edges from crossing between zones using different materials.

We optimize with two different approaches, one combinatorial and the other geometric. The combinatorial optimizer chooses which edges T uses in a fixed graph obtained from E CVT , E Len . It optimizes for E Ali , E Ori , E Zon . In this optimizer, the multiobjective equation parameterized by the cycle T is:

E Comb (T ) = E Zon + 𝛾 Comb (E Ali + E Ori ) (1)
The geometric optimizer takes this output as initialization and further optimizes T for E CVT , E Len , E Lap , E Ali , E Ori , while trying to keep E Zon unchanged. It directly manipulates the vertices and segments of T . The corresponding multi-objective equation is:

E Geo (T ) = E CVT + 𝛾 Lap E Lap + 𝛾 Len E Len + 𝛾 Obj (E Ali + E Ori ) (2)
Each optimizer performs the following optimization:

min T E (T ) (3)
with E being either E Comb or E Geo .

A detailed definition of each optimization term and its weights is provided in the Appendix, where we also derive the objective gradients, as these are required for the geometric optimizer.

There are several original aspects in the way we formulate our objectives. First, we formulate E CVT on the segments of T . This results in robust and precise optimization, naturally avoiding selfintersections without requiring a point resampling of T at every solver iteration. Equipped with the gradient, we optimize E CVT with the quasi Newtonian solver L-BFGS [LN89]. Second, we formulate the orientation objective in terms of the distribution of angles along T . By minimizing the distance to a uniform distribution, we obtain isotropic trajectories; by maximizing it, we encourage anisotropy -without choosing an a priori orientation. This leads to an optimal transport problem in the geometric solver context through the Wasserstein distance [START_REF] Villani | Topics in optimal transportation[END_REF]. Third, we consider the alignment objective by discrete integration of the target directions over the cells of the diagram of segments. This is achieved precisely and efficiently through a pre-computed table.

The objectives maintain their formulations in both the combinatorial and geometric optimizers.

Covering the surface with a graph

The first step of our approach is to construct a graph on which the combinatorial cycle optimizer will later operate. We seek to produce a planar graph embedded in P. We know from [GE04] that a Hamiltonian cycle can be efficiently computed in a bridgeless graph of degree three, where bridgeless means that the graph does not have any edge whose deletion would increase the number of connected components. Besides, we seek to obtain a graph with edges distributed evenly over P.

From these two requirements, we choose to create the graph from the vertices and edges of a Centroidal Voronoi Tessellation (CVT). We generate initial seeds by intersecting a triangular tiling with our polygon P (see Figure 3).

The spacing 𝑑 between adjacent vertices parameterizes the tiling. Each seed position is perturbed by a random offset, sampled in a disc of radius 𝑑 20 , in order to get rid of co-circular cases. The value of 20 was determined empirically. Tiling vertices inside P will be the initial seeds of our Voronoi diagram.

The spacing 𝑑 is the driving parameter for the curve interspace in the final result: with a lower 𝑑, the number of seeds 𝑆 increases, producing a longer cycle, which leads to a smaller curve interspace. However, choosing 𝑑 to obtain a specific interspace is not simple, as this also depends on the polygon P and the objectives we are minimizing (see Figure 13). In practice, 𝑑 is adjusted empirically by restarting the process with a lower 𝑑 when the interspace is too high and higher 𝑑 when the interspace is too low. After the Voronoi tessellation, the graph is composed of the Voronoi edges of the diagram, clipped by the polygonal contour whose boundary is also part of the graph, see Figure 4. A CVT graph is extremely likely to have nodes of degree three everywhere, as co-circular cases are extremely rare. Actually, we have never observed a node with a degree greater than three. The graph is bridge-less by construction, as each edge belongs to a cyclic cell. However, directly using a point CVT graph would yield poor results under alignment objectives. Indeed, during the combinatorial optimization, we select which edges belong to the cycle. If the user controls alignments through the vector field 𝑈, we can only hope to find a cycle with aligned edges if they exist in the graph. Therefore, the graph has to ideally contain some edges aligned along the non-zero vectors of 𝑈.

To align the edges, we introduce anisotropy in the Voronoi diagram. Instead of casting this as a metric distortion problem [LB13], submitted to COMPUTER GRAPHICS Forum (4/2022). we consider a diagram of many small segments -called needles.

P • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 𝑑 𝑑
The lengths and orientations of the needles induce cells that are elongated and oriented with respect to the control field 𝑈, see Figure 4. Each needle is centered on a supporting seed. The needle lengths and orientations are implicitly defined by the (point) Voronoi diagram of their supporting seeds. Let us consider the needle centered on a seed 𝑝, whose Voronoi cell we denote as C 𝑝 .

In order to obtain cell edges aligned with the vector field 𝑈, we want the needle to be aligned with the average of the direction field 𝑈 inside its own cell C 𝑝 . However, such a field has symmetry: a vector and its opposite represent the same direction. To break the symmetry, one can double the angle of 𝑈 (denoted by ∠𝑈), leading to a new vector field 𝑉. In that way, two opposite vectors in 𝑈 become equal in 𝑉, and integration inside a cell is then possible.

𝑉 = cos (2∠𝑈) sin (2∠𝑈) if 𝑈 ≠ (0, 0) 0 0 otherwise (4)
Now the angle of the needle can be expressed as 1 2 ∠V 𝑝 , half the angle of V 𝑝 , the integral of 𝑉 inside C 𝑝 , as defined in the appendix equation (32).

The needle length is computed as

L 𝑝 = | |V 𝑝 | | area( C 𝑝 ) • √ area( P)
𝜋𝑆 , where the square root term is a global scaling factor based on the total area covered with 𝑆 seeds and area(.) computes the area of a polygon. To ensure that needles never cross L 𝑝 is clamped so that the needle is fully contained in C 𝑝 .

We optimize for the seeds' positions under a segment CVT objective, using the L-BFGS algorithm [LN89]. We use the simplifying assumption that V 𝑝 undergoes only small changes at each iteration. Therefore, we only optimize seed positions, recomputing V 𝑝 , L 𝑝 at every iteration of the L-BFGS algorithm. The gradient of the segment CVT objective is computed at the needle extremities and used to compute the supporting seed positions' gradient. After convergence, we extract the graph as the edges of the segment Voronoi diagram and clip it with the boundary of P. The above procedure gives us an initial graph satisfying our requirements (see Figure 4).

Combinatorial optimization

Once a bridge-less planar graph with nodes of degree three is obtained, we proceed with the combinatorial cycle optimization. Our combinatorial optimizer's core principle is to start from a valid Hcycle and perform local operations on its edges, modifying the cycle while preserving its Hamiltonian property. These local changes allow us to perform an efficient stochastic exploration, optimizing the alignment, orientation, and zoning objectives. The use of local operations is inspired by recent graph theory results that reconfigure H-cycles in regular grids [START_REF] Nishat R | Bend complexity and Hamiltonian cycles in grid graphs[END_REF][START_REF] Nishat | Reconfiguration of Hamiltonian cycles and paths in grid graphs[END_REF]. We introduce a different set of operations as our graphs are irregular and use them in a stochastic optimization framework (Section 3.3.3).

The optimizer proceeds in three main steps. The first obtains a set of cycles fully covering the input graph but not yet connected in a global H-cycle (Section 3.3.1). For this, we rely on existing techniques. The second operates a graph refinement and reconnects the cycles into a single Hamiltonian cycle using local operations (Section 3.3.2). The graph refinement is a crucial step that enriches the local operations performed on cycles in general and H-cycles in particular. The third is the optimization process (Section 3.3.3).

Finding initial cycles

As a starting point, we follow the principles established in Gopi and Eppstein [START_REF] Gopi | Single-strip triangulation of manifolds with arbitrary topology[END_REF]. We compute a perfect matching on the graph [START_REF] Edmonds | Paths, trees, and flowers[END_REF]: a selection of edges such that every node is adjacent to exactly one selected edge. This can be done in linear time on bridge-less planar graphs of degree three [START_REF] Biedl | Efficient algorithms for Petersen's matching theorem[END_REF]. Then, a set of cycles covering all graph nodes is obtained by inverting this edge selection. However, these cycles are not connected in a single Hamiltonian cycle (see Figure 5). [START_REF] Gopi | Single-strip triangulation of manifolds with arbitrary topology[END_REF] reconnect the cycles through two strategies. The first is to identify specific configurations of edges allowing to split and reconnect neighboring cycles -however, such configurations are relatively rare in three-connected planar graphs. The second is to create double edges connecting cycles, which is always possible but locally creates a strong deviation from the even spatial distribution of edges we seek to achieve.

Gopi and Eppstein

Graph refinement and initial H-cycle

Our combinatorial optimizer relies on local modifications of a Hamiltonian cycle. We follow a similar approach to reconnect the initial cycles into a single H-cycle. The three-connected planar graph's low connectivity and mostly hexagonal faces make it challenging to perform simple modifications: any change impacts many edges at once. Therefore, we propose refining the graph, increasing the density of edges, and making it easier to find opportunities for local, simpler modifications. In the following, we refer to the faces of the three-planar graph as cells as they originated from Voronoi cells. We insert one additional node at the barycenter of every cell, splitting it into triangular faces. This offers a much richer set of opportunities to perform local manipulations of cycles.

After introducing a center node and splitting the cells into triangular faces in the graph (Figure 6), we update the cycles obtained from perfect matching to end up in a similar situation: a set of cycles covering all nodes.

We examine each cell in turn. In each cell, we choose one of the boundary edges (𝑎, 𝑏) that is part of the cycle. This edge is then split to go through the added center node 𝑛 as (𝑎, 𝑛) → (𝑛, 𝑏). The cycle now goes through 𝑛 in addition to (𝑎, 𝑏); see dashed lines in Figure 6, top right.

Note that all cells have at least two cycle edges along their boundary by construction. Indeed, we started from a perfect matching (Section 3.3.1) and inverted the selection, ensuring that every node is connected to two-cycle edges and a third edge not in the cycle. Thus, we cannot have two successive edges that do not belong to the cycle. Consequently, a cell made up of 𝑘 edges has at least 𝑘/2 of them in a cycle. Since a cell is made up of at least three edges, this implies that all cells have at least two edges in a cycle.

As we connect boundary cycle edges to cell centers, we may end up with a cell that no longer has any cycle edge along its boundary, making the new node 𝑛 unreachable (Figure 6, top right). When this happens, the cell steals this edge back from its neighbor, and the cell is tagged as frozen: no other cell can steal this edge from it. The neighbor now needs to connect its center again. It will either use another available boundary cycle edge or will have to steal it from another neighbor. The process propagates until all cases are resolved. This is illustrated in Figure 6. In all our experiments, the process terminates successfully, with minimal propagation. and 𝑒 11 and reconnect their endpoints with 𝑛 1 , 𝑛 3 and 𝑛 4 , respectively. However, we need to find an alternative since the resulting cycles do not cover the node 𝑛 2 . Bottom: here we remove edges 𝑒 5 and 𝑒 6 and reconnect their endpoints with nodes 𝑛 1 and 𝑛 2 respectively. This leads to a satisfactory result covering all nodes. Finally, once the new graph and cycles are obtained, we search for split-and-reconnect configurations to merge the cycles, as illustrated in Figure 7: an edge from a different cycle across two free edges. The much more densely connected graph provides a larger number of reconnection opportunities. The algorithm visits all cycle edges, searching for the candidate configuration. As soon as a split-and-reconnect operation is discovered, it is applied, joining the two cycles. The algorithm continues until no further split-andreconnect operation is possible.

At this stage, if multiple cycles remain, double edges are added, as illustrated in Figure 8. The algorithm goes along the cycles again, searching for the first free edge connecting another cycle and doubling it. Typically, only very few double edges are necessary, as experimentally verified in Table 1. In practice, doubling an edge is not so problematic as the geometric optimizer will quickly recover the proper interspace.

At this point, we have a valid initial H-cycle. This entire process is efficient: perfect matching is linear in the number of nodes, the cycle reconnection linear in the number of edges.

Stochastic H-cycle optimization

Objectives. Given an H-cycle, we evaluate the objectives E Ali , E Ori , and E Zon along the edge that belongs to the cycle. Please refer to Appendix A for their exact computations. However, some differences have to be clarified. For the alignment objective, as it would be very costly to compute a Voronoi diagram each time we update the score, the vector field integrals are no longer computed inside segment Voronoi cells but inside small rectangles surrounding edges. For the orientation objective, we track the angle distribution in a histogram having 𝑛 bins (𝑛 = 100 in our implementation). This histogram enables the computation of the objective value in constant complexity. The combinatorial optimization does not consider the surface-filling objectives (E CVT , E Len , E Lap ). The vertices are fixed, and the surface filling property is enforced by construction of the initial graph and the Hamiltonian nature of the cycle.

We update all objectives' current values as we locally modify the cycle by removing/adding edges. The objectives are weighted and combined as discussed in Section 3.1.

Evolving an H-cycle: local operations. The local operations we perform on the H-cycle are shown in Figure 9. There is a starting configuration and an updated configuration in each operation template. The starting configuration considers a set of nodes with a specific configuration of edges. The thick green edges are those that belong to the cycle; the thin black edges are those that do not currently belong to the cycle. Other edges may exist between the nodes and may belong to the cycle or not -these do not influence the operation.

Given a configuration that matches the template starting configuration, we can, in each case, modify the edges into the updated configuration while preserving the H-cycle. Each operation modifies the set of edges belonging to the cycle and thus the value of the objectives. We search for configurations by selecting a node at random and then considering the edges surrounding it. Matching the templates can be achieved in constant time for all operations but one (rightmost in Figure 9). A global connectivity check is required between two of the nodes, indicated in Figure 9 by a dashed green curve. We perform it in constant time using an auxiliary data structure. However, this data structure requires a linear time update after each actual modification of the graph -these only happen when a local change is accepted, typically when it improves the objective function's value.

Optimization loop. Our algorithm follows a standard stochastic optimization framework. Starting from an initial H-cycle, we produce a population of 𝐾 × 𝐶 cycles, cloning the initial one. Each of these H-cycles is evolved with the shuffle algorithm 1.

Algorithm 1 The shuffle algorithm applies local operators to evolve and explore input H-cycle improvements with respect to our objectives.

algorithm Shuffle(HCycle T) { for (shuffle = 0 to ShuffleIter) {

T = applyRandomOp(T) if (E Comb ( T ) <= E Comb ( T)) { T = T } else if (E Comb ( T ) < E Comb ( T) + ThreshHigh)) { if ((1 + 𝑁 * (E Comb ( T ) -E Comb ( T))) * rand() < ThreshLow) T = T } else if (shuffle % Mod == 0) { T = T } } } for (improve = 0 to ImproveIter) { T = applyRandomOp(T) if (E Comb ( T ) <= E Comb ( T)) { T = T } } return T; }
We next select 𝐾 champions from the 𝐾 × 𝐶 cycles, promoting diversity. We measure the distance between two H-cycles as the number of edges they do not have in common. We initialize the set of champions with the H-cycle having the best objective. Then we iteratively add the H-cycles that are the most distant from the already selected champions (Hochbaum-Shmoys algorithm [START_REF] Hochbaum | A best possible heuristic for the k-center problem[END_REF]) until 𝐾 are selected.

We then enter the main optimization loop, where champions are iteratively improved. At each iteration, each champion is cloned into 𝐶 candidates evolved with Algorithm 1. The champion is replaced by the best of its candidates, but only if it improves the objective compared to the champion. This is performed for a fixed number of iterations 𝐼. In practice, we use ShuffleIter=ImproveIter=8𝑁, Thresh-High=0.5, ThreshLow=1, Mod=400 and 𝐶 = 2, 𝐾 = 12, 𝐼 = 3. Here, 𝑁 is the number of nodes inside the cycle. The algorithm uses some randomness through the function rand() that returns a uniform number between 0 and 1. Finally, this is easily parallelized; we use one thread per champion.

Cycle geometry optimization

The input of the geometry optimization is the H-cycle T obtained by the combinatorial optimizer. Our goal is now to optimize the position of the vertices of T in order to obtain smooth trajectories evenly distributed while following all our other objectives (alignment, orientation, zoning).

Initialization An important question for the geometric optimizer is which length to target for the cycle with E Len . Let us recall that the total length indirectly controls the final spacing between folds of the trajectory.

The cycle produced by the combinatorial optimizer has a global directionality in agreement with the alignment and orientation objectives. However, it also has many local corners (high frequencies) coming from the initial tessellation. Now that the cycle connectivity is determined, the geometric optimizer can smooth out these higher frequencies. The smoothing objective E Lap achieves this. However, the high frequencies coming out of the combinatorial optimizer artificially lengthens the cycle. Therefore, we apply a few explicit smoothing steps before computing the target length. This results in a filtered cycle whose length is more representative of our spacing objective.

We apply 5 times the following smoothing operator to each vertex 𝑣 𝑣 𝑣 𝑖 , while ensuring that no self-intersections appear:

𝑣 𝑣 𝑣 𝑖 ← 𝑣 𝑣 𝑣 𝑖 + 1 10 𝑣 𝑣 𝑣 𝑖 + + 𝑣 𝑣 𝑣 𝑖 --2𝑣 𝑣 𝑣 𝑖 (5)
After applying (5), the target length 𝐿 0 becomes the length of T . The geometric optimizer starts from this smoothed cycle.

Optimization We then minimize the geometric objective E Geo defined in [START_REF][END_REF]. We encourage the reader to have a look at the different terms of this objective in the appendix. Our geometric optimizer is supported by a quasi Newtonian solver L-BFGS [LN89]. This is an Before smoothing (5) After smoothing (5)

After optimization iterative solver that only needs to compute the value of the objective and the gradient or at least an approximation of it to minimize the objective. The solver then approximates the Hessian using the gradient of the last 𝑀 iterations. 𝑀 is fixed before optimization and has a value generally less than 10; we use 7 in our implementation.

Figure 10 illustrates the impact of the initial smoothing (5) and the optimizer over the cycle geometry.

Computing the per-segment area After the process has converged, we are ready to compute the per-segment area along T . Once multiplied by a thickness, it represents the volume of material that has to be extruded around the segment. This enables the deposition to adapt to each segment's surroundings (see Section 4.2). The area is computed from the segment's Voronoi cell, which is added half the area of the point Voronoi cells of its extremities.

Results

This section presents numerical evaluations and experimental results, including plates 3D-printed with multiple materials.

Numerical evaluations

Figure 11 shows the evolution of the two optimization stages (combinatorial and geometrical). As can be seen, both optimization solvers operate effectively to reduce the value of the objectives.

For the geometrical optimization, 𝛾 Len E Len acts as a regularizer of the total length of T , and the optimizer attempts to maintain its value. Figure 12 shows the effect of disabling the optimization of global terms (E Ori , E Ali ). As can be seen, these terms are critical in obtaining a trajectory well aligned with the target field. Step Weighted objective Figure 13 reveals the link between the parameter 𝑑 controlling the initial spacing between adjacent seeds (see Section 3.2), and the interspace 𝜏 encountered along the cycle. The interspace 𝜏 is computed from the Voronoi cells, taking the distance between both sides across the trajectory. The relationship between 𝑑 and the final interspace behaves consistently across examples and can be predicted reasonably well using the curve of Figure 13. We can also observe that the variation across models is slight (within ±0.1 mm of each other's).

E Geo E CVT 𝛾 Lap E Lap 10 • 𝛾 Len E Len 50 • 𝛾 Obj E Ali 50 • 𝛾 Obj E Ori
In addition, Figure 14 shows that the interspace values follow a (tilted) Gaussian distribution. This simplifies pushing a varying flow (by adjusting motion speed) during fabrication. Interestingly we observe a change in statistics between the objective. However, they all remain within small bounds, as their means are all within 0.05 mm of each other. One can also notice that the variance of the anisotropy is lower, as it is easier to obtain a constant interspace with zigzags.

Figure 15 demonstrates the effect of the Laplacian term 𝛾 Lap (see Section A). While both results are compelling, using the Laplacian term avoids sharp turns, effectively smoothing the cycle.

Figure 16 shows some timing results. Both optimizers can execute iterations in a reasonable time, with a roughly polynomial correlation with respect to 𝑑.

3D printed results

We carried out 3D printing of parts optimized with our method. We employed filament fused fabrication with PLA.

Multi-material illustrations

We fabricated multi-material plates with complex trajectory patterns. These prints were fabricated using a constant deposition width to highlight the trajectory better. We typically print five layers. Figure 16: Timings of the brain model (see Figure 11). Left: computing time of the function Shuffle (Fig. 1) with respect to 𝑑. Right: computing time of a single iteration of the L-BFGS solver in the geometric step.

We use a 3-way Diamond nozzle by Reprap.me. This nozzle takes three different filaments as input, mixed into the nozzle chamber, and exit from a single hole. The mixture is passive, so the materials are co-extruded in a toothpaste-like pattern at a small scale. Our printer is a Prusa-like assembly, with three feeders for each of the three filaments. The printer firmware accepts special GCode instructions defining the mixing ratio at each extrusion move. In our experience transitioning from one mixture to another requires pushing at least 6.1𝑚𝑚 3 , resulting in a relatively long transition, e.g., 76 mm of deposition trajectory for a 0.4 mm nozzle at 0.2 layer height. In this context, the zoning and alignment objectives are important to minimize the number of transitions and align them along zone boundaries. We printed at a speed of 10 mm/sec, using colored translucent PLA filament.

Figure 17 shows the brain model 3D printed with a variety of materials (colors) under different objectives combined (see Figure caption for details). Two versions were printed with different spacings 𝑑. As can be seen, the 3D printed versions closely follow the input specifications, with clearly defined zones both in material and alignment/orientation objectives. The single continuous deposition produces an even, reliable material flow, leading to nearly optimal operational conditions for the extrusion process.

The leaf model shown printed in Figure 18 has many narrow regions, and an internal skeletal structure. This shows that despite being a single cycle, the trajectory can fill narrow regions orthogonally to their main directionality -this would be difficult to achieve with growth techniques that have to fit an even number of trajectories everywhere (see Section 2). The skeleton is well captured despite the multi-material transitions.

Figures 19 and20 show other examples of shapes with narrow spaces and multi-material transitions. The batman logo uses an alignment objective parallel to the boundaries. The zebra uses anisotropy for the background and alignment along boundaries within the zebra stripes.

Dense plates

We compared two methods on an isotropic cycle filling a square to evaluate how dense our prints can be. Isotropy results in a high number of turns. The first method uses our persegment area computation, varying the amount of deposited material along the cycle. Note that the material's feed rate was kept constant (the speed of the motor pushing the filament) while the XY motion speed was varied. The second method extrudes at a constant speed an average amount of material, computed from the plate's total area -this is equivalent to considering that the average 'width' along the cycle is the total area divided by cycle length. The print speed is slightly impacted by the type of infill due to the acceleration limits. Comparing on a CR10S, printing a 40x40 mm square with a speed target of 20mm/sec, the isotropic infill reaches an average deposition speed of 19.75 mm/sec versus 19.94 mm/sec for an anistropic one.

As shown in Figure 21 the per-segment area achieves a denser result, leaving a few gaps. These examples are fabricated at 0.2 mm layer height, with two layers, using a 0.6 mm nozzle at an average speed of 20 mm/sec on a Prusa i3 printer. They are printed in sequence and at the same bed location. Both use the same total volume of material: 281 𝑚𝑚 3 for the squares, 1955 𝑚𝑚 3 for the leaves. The average spacing around the cycle is around 1.3 mm. Note that we printed at an increased temperature (215 degrees) to allow the molten plastic to spread more easily. The nozzle width must also be within reasonable margins of the nominal spacing (obtained as total area divided by cycle length). 3D parts Beyond plates, we print 3D objects using our technique. Figure 22 compares three printed versions of the CuteOcto model, where our method is used in isolation and compared with a standard zigzag infill, also used in isolation. We print the CuteOcto first using our isotropic fill and next using a vector field that encourages the deposition to conform to the surface. In both cases, we obtain a solid object that prints with no travel moves beyond those strictly necessary (layers with two disjoint areas, e.g., front paw).

As we use the infill in isolation, we can see that the turns along the surface produce small depressions. However, even in the isotropic case, the surface quality remains better than that produced by a standard zigzag infill, creating turns everywhere along the boundary. Using the vector field (Figure 22, middle) further encourages the paths to circulate along the part boundary.

In standard practice, to obtain a smooth surface, an external contour (perimeter) is used, with the infill hidden behind. Since our infill is a cycle, and since external contours would also be cyclic, it would be simple to reconnect everything in a global cyclic trajectory, still producing a continuous extrusion.

Our technique is also useful on multi-material parts, as shown in Figure 23. Here, we use zoning to transition a cyclic isotropic infill between two material mixtures. The extrusion is continuous, and the transition between mixtures occurs in a user-specified zone (Figure 23, bottom left). In contrast, a two nozzle print on an Ultimaker 3 requires switching every layer and printing an extra tower to prime the nozzles, wasting material and time. Finally, an advantage of using mixtures is to print with none of the original colors: as long as the mixture is stable (e.g., 30% -70% mix), a consistent color is produced [SMB * 19].

Performance

We measure performance on a Core i7-10750H laptop (2.6GHz) with 16 GB of RAM. Our implementation is a non-optimized research prototype. We report in Table 2 some timings. The combinatorial optimizer runs several threads in parallel (12). On CuteOcto, timings are the total for all layers (noting that CuteOcto has a large base followed by smaller layers above). Init and geometric optimizers are parallelized when multiple slices are computed. In all cases, the total computation time is below fabrication time. The algorithm could be streamed as slices are processed and sent to the printer and multiple slices processed on different cores. The stopping criteria offer a time-quality tradeoff (see Figure 16). We use two settings: one for cases with zoning and one for single objective cases (the combinatorial optimizer typically converges faster on these).

In terms of complexity, initializing the graph is 𝑂 (𝑁), where 𝑁 is the number of cycle vertices, 𝑁 = 𝑂 (

𝑎𝑟 𝑒𝑎 ( 𝑃) 𝑑 2
). Each iteration of the combinatorial optimizer is 𝑂 (𝑁 2 ). The 𝑁 2 is due to the worst linear case of updating the connectivity check data structure of the transpose operator (Figure 9, rightmost). Each iteration of the geometric optimizer iteration is 𝑂 (𝑁𝑙𝑜𝑔(𝑁)). Table 2: Timings of different models (averaged, five runs).

Discussion

In this section, we provide comparisons and additional discussion. Fermat spiral. The method of [ZGH * 16] constructs a concentric spiral explicitly. Our approach produces a trajectory that is well aligned with the concentric field and also exhibits a more even inter-spacing. In particular, the offsetting used during the Fermat spiral construction tends to accumulate gaps along the medial axis, a defect that would repeat and align throughout layers. Thanks to their randomized nature, our trajectories are less subject to defect accumulation across layers. We can also produce a cyclic trajectory following the direction orthogonal to the concentric field (the green region in Figure 25, see also Figure 12 left), which cannot be achieved with Fermat spirals.

Comparisons

Labyrinths and mazes. Our technique produces an outcome similar to that of Pedersen and Singh [PS06], following the field closely. It can also be seen that our result crosses between the three zones four times (two between green and orange -once around the arms, once around the feet -and two between orange and gray) which is the minimum possible. A growth approach will always cross at least six times (see also Figure ??).

Limitations of growth approaches

Growth approaches offer multiple controls, some of which are similar to our goals, anisotropy, and zoning in particular. However, we are not aware of any method explicitly targeting isotropy.

In particular, the method of Pedersen and Singh [PS06] demonstrates impressive artistic results using these controls. This technique is primarily designed as an interactive authoring tool. Some of its properties do not translate to the automatic generation of trajectories in a large number of layers, as we discuss now. Figure 22: CuteOcto model 3D printed with our isotropic (left), our vector concentric fills (middle), and standard zigzag for comparison (right). Please note that in standard prints, the infill 'sides' are usually hidden beneath an external contour parallel perimeter, but in this example we compared the infills in isolation.

In the discussion below, we compare the number of iterations between methods. The growth process starts from a circle, and our method starts from our combinatorial initialization that takes 𝑂 (𝑁) with 𝑁 number of points along the trajectory. Our geometric optimizer and the growth process of Pedersen and Singh both have an 𝑂 (𝑁 log(𝑁)) complexity per iteration. However, our approach performs a quasi-Newton optimization that converges towards a (local) minima in a few iterations. In contrast, the growth process does not converge: in our experience, the result keeps moving even after invading the entire domain. For the growth process, we used the implementation https://github.com/twentylemon/ organic-labyrinth that we improved and extended.

A drawback of growth processes is to take significant amounts of time to evolve the curve, especially as the growth rate has to be limited to a fraction of the target inter-spacing to guarantee a stable growth [PS06]. To alleviate this issue, the growth is first performed with a larger inter-space that is then progressively reduced, as shown in Figure 26a. This is done manually in an interactive setting. In Figure 26a we decreased the interspace at manually selected points. In contrast, our technique immediately starts at the chosen inter-space. As a result, in Figure 26b our approach converges in 170 iterations versus 1700 for the multi-scale growth process. This advantage is further increased in geometries with narrow passages. These hinder the use of a decreasing inter-space: with a large inter-space, the curve cannot invade past the narrow regions. This is illustrated in Figure 26c. Here the growth has to start with an inter-space smaller than the minimal feature size to reach every location. A more elaborate initialization -beyond a simple circle -would alleviate this issue. However, doing so in a general setting where the shape contains holes is not trivial and is what our combinatorial optimizer provides. As a result, in Figure 26d our approach converges in 260 iterations versus 110000 for the growth process. Please also observe that continuously growing from a simple closed curve limits the cycles that can be achieved. For instance, this implies there will always be an even number of trajectories between two outer contours, making trajectories such as the closeup shown in Figure 18 unfeasible. This further limits following specific orientations in narrow features, such as the folds in Figure 18 which would not be possible with two trajectories side by side.

Finally, it is worth noting that after combinatorial initialization, we could use the optimizer of Pedersen and Singh [PS06]. However, the dynamic evolution process of the curve relies on a number of parameters that are difficult to tune to achieve specific properties (e.g., anisotropy, isotropy). This is not a major issue in an interactive setting where adjustments are performed live but proves difficult for trajectory computations where the parameters cannot be easily adjusted for every single situation. In contrast, our curves are globally optimized.

Influence of orientation on mechanical properties

Multiple factors affect the mechanical properties of the fabricated layers. Among the most important are the bonding strength, extrusion temperature, printing speed, extrusion rate, nozzle and filament diameter, viscosity, and thermal expansion coefficient of the deposition material [ATB * 18, CdSC19]. A tightly controlled manufacturing environment and protocols are necessary to obtain consistent mechanical results.

Fortunately, several prior works studied this and indicate that the orientation of the deposition paths plays a role in the mechanical properties of each layer. For instance, Villacres et al. [START_REF] Villacres | Additive manufacturing of shape memory polymers: effects of print orientation and infill percentage on mechanical properties[END_REF] found the orientation of parallel paths with an infill density close to 100% has a monotonic relationship with stiffness: Young's modulus is highest for the parallel direction and lowest for the orthogonal direction. Recent studies have also been interested in producing paths that exhibit transversely isotropic mechanical properties, as they offer mechanical advantages compared to the traditional parallel paths [LXM * 19]. Other examples of works using the mechanical influence of trajectory orientations are [SIM16, TM17, LY17, KWW19a, XLM20, FZZ * 20].

A promising direction is to increase the influence of orientation, for instance, using on-purpose non-optimal deposition flow and temperature. This would potentially make bonding between neighboring paths weaker, increasing flexibility in these directions.

Limitations

To obtain denser plates, our approach requires printing while controlling motion speed for a varying flow deposition (as has been successfully demonstrated in recent works [HKD * 20]). However, some tiny porosities remain in challenging cases (see the leaf in Figure 21).

Our current graph construction tends to create edges that take many orientations, even when biased by the needles. This puts more pressure on the geometric optimizer. Other graph constructions, perhaps inspired by quad-meshing, could alleviate this issue.

Conclusions

Our technique offers control over material deposition orientation while manufacturing a solid layer along a strictly continuous extrusion trajectory. It globally optimizes a space-filling trajectory for isotropy, anisotropy, vector field alignment, and zoning. As we have demonstrated, this is of particular interest for multi-material prints and aesthetic purposes. Continuity is also an essential consideration for clay, ceramics, silicon, and concrete materials.

Our framework also supports a wide range of optimization objectives. We believe this will prove helpful for Additive Manufacturing processes using challenging materials, such as the fabrication of plates having isotropic elastic behavior [LXM * 19]. To obtain the gradient with respect to the initial path's vertices, we will have to make the following chain rule update at the end of our computations: We consider computing the objective and its gradient inside a single zone in the following. For simplicity and not introducing new notations, we suppose the zone to be P, and we still use the notation T for the intersection of the path inside the zone. Notice that T could no longer be a cycle but a set of paths having two endpoints. In the same way, we still denote by 𝑣 𝑣 𝑣 𝑖 the path's vertices, including the ones that intersect the boundary of the zone.

G (𝑖) CVT ← G (𝑖) CVT + (1 -𝑡) • G (𝑢) CVT G ( 𝑗) CVT ← G (𝑖) CVT + 𝑡 • G
The segment Voronoi diagram over T gives the following partition of P into two types of Voronoi cells:

• Point cell A cell C 𝑖 for each vertex 𝑣 𝑣 𝑣 𝑖 of T , containing all points of P that are nearest to 𝑣 𝑣 𝑣 𝑖 than to any other point of T :

C 𝑖 = {𝑥 ∈ P | ∀𝑦 ∈ T , 𝑥 -𝑣 𝑣 𝑣 𝑖 𝑥 -𝑦 } (11)
• Segment cell A cell C 𝑖 𝑗 for each open segment ]𝑣 𝑣 𝑣 𝑖 , 𝑣 𝑣 𝑣 𝑗 [ of T that contains all points of the shape P that are nearest to this segment than to any other point of T :

C 𝑖 𝑗 = 𝑥 ∈ P | ∀𝑦 ∈ T , dist 𝑥, ]𝑣 𝑣 𝑣 𝑖 , 𝑣 𝑣 𝑣 𝑗 [ 𝑥 -𝑦 (12)
This partition is illustrated Figure 27. Once this partition is obtained, the distance to the path can be obtained as a close form inside each cell. The CVT objective is rewritten as:

E CVT = 𝑁 ∑ 𝑖=1 ∫ C 𝑖 𝑥 -𝑣 𝑣 𝑣 𝑖 2 𝑑𝑥 + ∑ (𝑖 𝑗) ∈T 𝐸 ∫ C 𝑖 𝑗 dist 𝑥, ]𝑣 𝑣 𝑣 𝑖 , 𝑣 𝑣 𝑣 𝑗 [ 2 𝑑𝑥 = 𝑁 ∑ 𝑖=1 E CVT,𝑖 + ∑ (𝑖 𝑗) ∈T 𝐸 E CVT, (𝑖 𝑗)
(13) Where T 𝐸 = {(𝑖, 𝑖 + ) | 𝑖 = 1, ..., 𝑁 } is the set of edges of the cycle. The boundary of a Voronoi cell consists of bisectors sections between two sites (which can be a vertex or an open segment). All the bisectors are straight lines apart from the parabolic bisectors between a vertex and an open segment. In order to simplify computations, we approximate the section of a parabola with 𝑘 straight line segments. In our implementation, we choose 𝑘 = 10, which provides a fine tessellation.

We present how we derive the objective and gradient terms in the following. From now on, C 𝑖 and C 𝑖 𝑗 are assumed to be composed only of straight-line segments due to the parabolic bisector sections' approximation.

submitted to COMPUTER GRAPHICS Forum (4/2022).

Point cell

We want to derive the objective E CVT,𝑖 associated with one point cell C 𝑖 . First, we rewrite it as:

E CVT,𝑖 = ∫ C 𝑖 𝑥 -𝑣 𝑣 𝑣 𝑖 2 𝑑𝑥 = ∫ C 𝑖 𝑥 2 𝑑𝑥 -2 ∫ C 𝑖 (𝑥.𝑣 𝑣 𝑣 𝑖 )𝑑𝑥 + ∫ C 𝑖 𝑣 𝑣 𝑣 𝑖 2 𝑑𝑥 = ∫ C 𝑖 𝑥 2 𝑥 + 𝑥 2 𝑦 𝑑𝑥 -2 ∫ C 𝑖 𝑣 𝑣 𝑣 𝑖, 𝑥 𝑥 𝑥 + 𝑣 𝑣 𝑣 𝑖,𝑦 𝑥 𝑦 𝑑𝑥 + ∫ C 𝑖 𝑣 𝑣 𝑣 𝑖 2 𝑑𝑥
(14) where 𝑥 𝑥 and 𝑣 𝑣 𝑣 𝑖, 𝑥 are, respectively, the 𝑥 coordinate of the point 𝑥 and the cycle vertex 𝑣 𝑣 𝑣 𝑖 (and the same applies for the subindex 𝑦).

The gradient of E CVT,𝑖 with respect to 𝑣 𝑣 𝑣 𝑖 𝑖 𝑖 is:

G 𝐶𝑉 𝑇 ,𝑖 = 2 𝑣 𝑣 𝑣 𝑖 𝑖 𝑖 ∫ C 𝑖 1𝑑𝑥 - ∫ C 𝑖 𝑥𝑑𝑥 (15)
New terms that appear in Equations 14 and 15 are moments of order 𝑛 = 0, 1, 2 of the point cell. The Shoelace formula allows the computation of the signed area of a polygon (the moment of order 0) and generalizes to moments greater than 𝑛 > 0.

We denote by 𝑝 𝑗 the vertices of the polygon C 𝑖 and by ( 𝑗 𝑘) its edges oriented in clockwise order. For simplicity, we write ( 𝑗, 𝑘) ∈ C 𝑖 to iterate over the edges of the cell. Shoelace formula gives us:

∫ C 𝑖 𝑥 𝑛 𝑦 𝑑𝑥 = 1 (𝑛 + 1) (𝑛 + 2) ∑ ( 𝑗 𝑘) ∈ C 𝑖 ( 𝑝 𝑘, 𝑥 -𝑝 𝑗, 𝑥 ) 𝑛+1 ∑ 𝑙=0 𝑝 𝑙 𝑘,𝑦 𝑝 𝑛+1-𝑙 𝑗,𝑦 (16)
∫ C 𝑖 𝑥 𝑛 𝑥 𝑑𝑥 = 1 (𝑛 + 1) (𝑛 + 2) ∑ ( 𝑗 𝑘) ∈ C 𝑖 ( 𝑝 𝑗,𝑦 -𝑝 𝑘,𝑦 ) 𝑛+1 ∑ 𝑙=0 𝑝 𝑙 𝑘, 𝑥 𝑝 𝑛+1-𝑙 𝑗,𝑥 (17)
Segment cell For (𝑖, 𝑗) ∈ T 𝐸 , we want to compute the objective associated with the segment cell C 𝑖 𝑗 . First, we project our points on a new orthogonal basis to facilitate subsequent derivations. This basis is composed of the two following unit vectors v and ŵ:

v = 𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 , ŵ = 𝑅 𝜋/2 v (18)
Where 𝑅 𝜋/2 is the matrix of a rotation by an angle 𝜋/2. We then consider the following affine isometric transformation 𝑚:

𝑚 : 𝑥 ↦ → (𝑥 -𝑣 𝑣 𝑣 𝑖 ). v (𝑥 -𝑣 𝑣 𝑣 𝑖 ). ŵ (19)
Two functions characterize the transformed cell 𝑚 C 𝑖 𝑗 : 𝑓 < 0 and 𝑔 > 0 representing the lower border and the cell's upper border. They are illustrated Figure 28. The support used for these two functions is the interval ]0, 1[. We then rewrite the segment objective as:

E CVT, (𝑖 𝑗) = ∫ C 𝑖 𝑗 dist 𝑥, ]𝑣 𝑣 𝑣 𝑖 , 𝑣 𝑣 𝑣 𝑗 [ 2 𝑑𝑥 = 𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 ∫ 1 0 ∫ 𝑔 (𝑡) 𝑓 (𝑡) 𝑦 2 𝑑𝑦𝑑𝑡 (20)
Gradients with respect to the two endpoints 𝑣 𝑣 𝑣 𝑖 , 𝑣 𝑣 𝑣 𝑗 are given by: For each vertex 𝑝 𝑘 of the cell C 𝑖 𝑗 , we let 𝑞 𝑘 = 𝑚( 𝑝 𝑘 ) be the associated vertex in the new transformed cell. Using the expression of the 𝑦 moment of order 2 (see Equation 16), we get:

G (𝑖) CVT, (𝑖 𝑗) = -𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 ŵ ∫ 1 0 (1 -𝑡) 𝑔(𝑡) 2 -𝑓 (𝑡) 2 𝑑𝑡 (21) G ( 𝑗) CVT, (𝑖 𝑗) = -𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 ŵ ∫ 1 0 𝑡 𝑔(𝑡) 2 -𝑓 (𝑡) 2 𝑑𝑡 (22)
E CVT, (𝑖 𝑗) = 1 12 ∑ (𝑘𝑙) ∈ C 𝑖 𝑗 𝑞 𝑙,𝑥 -𝑞 𝑘,𝑥 3 ∑ 𝑚=0 𝑞 𝑚 𝑘,𝑦 𝑞 3-𝑚 𝑙,𝑦 (23)
For the gradient with respect to the two endpoints we first introduce new values:

𝑢 𝑘 = 𝑞 𝑘,𝑥 𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 , 𝑢 𝑙 = 𝑞 𝑙,𝑥 𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 (24)
and we then have:

G (𝑖) CVT, (𝑖 𝑗) = ŵ ∑ (𝑘𝑙) ∈ C 𝑖 𝑗 𝑞 𝑘,𝑥 -𝑞 𝑙,𝑥 12 3 ∑ 𝑚=1 (4 -𝑚𝑢 𝑘 + (4 -𝑚)𝑢 𝑙 ) 𝑞 𝑚-1 𝑘,𝑦 𝑞 3-𝑚 𝑙,𝑦 (25)
G ( 𝑗) CVT, (𝑖 𝑗) = ŵ ∑ (𝑘𝑙) ∈ C 𝑖 𝑗 𝑞 𝑘,𝑥 -𝑞 𝑙,𝑥 12 3 ∑ 𝑚=1 (𝑚𝑢 𝑘 + (4 -𝑚)𝑢 𝑙 ) 𝑞 𝑚-1 𝑘,𝑦 𝑞 3-𝑚 𝑙,𝑦 (26)

Length preservation objective

We denote by 𝐿 0 the target length. The length preservation objective keeps the length 𝐿 of T as close as possible to 𝐿 0 :

E Len = (𝐿 -𝐿 0 ) 2 = ∑ (𝑖 𝑗) ∈T 𝐸 𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 -𝐿 0 2 (27)
The gradient with respect to 𝑣 𝑣 𝑣 𝑖 is given by:

G (𝑖) Len = 2 (𝐿 -𝐿 0 ) 𝑣 𝑣 𝑣 𝑖 -𝑣 𝑣 𝑣 𝑖 + 𝑣 𝑣 𝑣 𝑖 -𝑣 𝑣 𝑣 𝑖 + + 𝑣 𝑣 𝑣 𝑖 -𝑣 𝑣 𝑣 𝑖 - 𝑣 𝑣 𝑣 𝑖 -𝑣 𝑣 𝑣 𝑖 - (28)

Smoothing objective

In order to have a smoother trajectory T we minimize the following Laplacian:

E Lap = 𝑁 ∑ 𝑖=1 𝑣 𝑣 𝑣 𝑖 - 1 2 𝑣 𝑣 𝑣 𝑖 + + 𝑣 𝑣 𝑣 𝑖 - 2 (29)
This objective smooths the sharp angles of T , redistributing its vertices. The gradient with respect to 𝑣 𝑣 𝑣 𝑖 is:

G (𝑖) Lap = 3𝑣 𝑣 𝑣 𝑖 + 1 2 𝑣 𝑣 𝑣 𝑖 ++ + 𝑣 𝑣 𝑣 𝑖 ---2 𝑣 𝑣 𝑣 𝑖 + + 𝑣 𝑣 𝑣 𝑖 - (30)
Normalized weights used for the two previous objectives are:

𝛾 Len = 0.4 • area(P)/𝑁 2 𝛾 Lap = 0.3 • area(P)/𝑁 (31)
submitted to COMPUTER GRAPHICS Forum (4/2022).

Alignment objective

Recall that we compute a vector field 𝑉 by doubling the angle of 𝑈 (Section 3.2). We then denote by V 𝑖 the integral of 𝑉 over C 𝑖 and by V 𝑖 𝑗 the integral of 𝑉 over C 𝑖 𝑗 :

V 𝑖 = ∫ C 𝑖 𝑉 (𝑥)𝑑𝑥 V 𝑖 𝑗 = ∫ C 𝑖 𝑗 𝑉 (𝑥)𝑑𝑥 (32)
The objective function E Ali seeks to minimize the angle difference between T and the integrated vector field inside each cell:

E Ali = 𝑁 ∑ 𝑖=1 V 𝑖 𝑍 ∠ 𝑣 𝑣 𝑣 𝑖 + -𝑣 𝑣 𝑣 𝑖 -- ∠V 𝑖 2 2 + ∑ (𝑖 𝑗) ∈ T 𝐸 V 𝑖 𝑗 𝑍 ∠ 𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 - ∠V 𝑖 𝑗 2 2 (33)
Angles can be defined using arctan function. We also write down its gradient for later:

∠𝑣 =        arctan(𝑣 𝑦 /𝑣 𝑥 ) if 𝑣 𝑥 > 0 arctan(𝑣 𝑦 /𝑣 𝑥 ) + 𝜋 if 𝑣 𝑥 < 0 (𝜋𝑣 𝑦 ) / (2|𝑣 𝑦 |) if 𝑣 𝑥 = 0 ∇∠𝑣 = 𝑅 𝜋/2 𝑣 𝑣 2 (34)
The values of ∠V 𝑖 and ∠V 𝑖 𝑗 in (33) are taken such that the absolute differences between angles are less than 𝜋/2. Finally, 𝑍 allows us to normalize the sum to obtain a weighted mean of squared angle differences. We have:

𝑍 = 𝑁 ∑ 𝑖=1 V 𝑖 + ∑ (𝑖 𝑗) ∈ T 𝐸 V 𝑖 𝑗 ∫ P 𝑉 (𝑥) 𝑑𝑥 -1 (35)
The integral of the norm over P is the total area in which the vector field is non-zero. This is useful for adequately weighing this objective with the isotropic/anisotropic one.

The gradient of this objective is approximate since the variations of the integrals V 𝑖 , and V 𝑖 𝑗 are ignored. The gradient with respect to 𝑣 𝑣 𝑣 𝑖 is approximated as follows:

G (𝑖) Ali = V 𝑖 - 𝑍 𝑅 𝜋/2 (𝑣 𝑣 𝑣 𝑖 -𝑣 𝑣 𝑣 𝑖 --) 𝑣 𝑣 𝑣 𝑖 -𝑣 𝑣 𝑣 𝑖 -- 2 ∠ 𝑣 𝑣 𝑣 𝑖 -𝑣 𝑣 𝑣 𝑖 --- 𝑖 - 2 - V 𝑖 + 𝑍 𝑅 𝜋/2 (𝑣 𝑣 𝑣 𝑖 ++ -𝑣 𝑣 𝑣 𝑖 ) 𝑣 𝑣 𝑣 𝑖 ++ -𝑣 𝑣 𝑣 𝑖 2 ∠ 𝑣 𝑣 𝑣 𝑖 ++ -𝑣 𝑣 𝑣 𝑖 - ∠V 𝑖 + 2 + V 𝑖 -𝑖 𝑍 𝑅 𝜋/2 (𝑣 𝑣 𝑣 𝑖 -𝑣 𝑣 𝑣 𝑖 -) 𝑣 𝑣 𝑣 𝑖 -𝑣 𝑣 𝑣 𝑖 - 2 ∠ 𝑣 𝑣 𝑣 𝑖 -𝑣 𝑣 𝑣 𝑖 -- ∠V 𝑖 -𝑖 2 - V 𝑖𝑖 + 𝑍 𝑅 𝜋/2 (𝑣 𝑣 𝑣 𝑖 + -𝑣 𝑣 𝑣 𝑖 ) 𝑣 𝑣 𝑣 𝑖 + -𝑣 𝑣 𝑣 𝑖 2 ∠ 𝑣 𝑣 𝑣 𝑖 + -𝑣 𝑣 𝑣 𝑖 - ∠V 𝑖𝑖 + 2 (36

Orientation objective

We control the probability distribution of the orientations along T to be either as isotropic or as anisotropic as possible. To do so, we minimize/maximize the squared 2-Wasserstein distance [START_REF] Villani | Topics in optimal transportation[END_REF] between a uniform distribution U over [0, 𝜋[ and the distribution of the orientations in T that we denote D T . The orientation distribution of T can be defined as: 𝐹 𝛼 is defined as 𝐹 -𝛼 and the pseudo inverses (𝐹 𝛼 ) (-1) and 𝐺 (-1) T are given by:

(𝐹 𝛼 ) (-1) (𝑡) = 𝜋(𝑡 +𝛼), 𝐺 In the distance (43), the infimum is obtained at the value 𝛼 0 defined as the mean of 1 𝜋 𝐺 (-1) T (𝑡) -𝑡 over [0, 1], which is:

𝛼 0 = 1 𝜋 E[D T ] - 1 2 (47)
We then obtain: To simplify the following expressions, we let:

Θ 𝑘 𝑖 𝑗 = 𝜋𝑇 𝑘 𝑖 𝑗 + 𝜋𝛼 0 (49)
Finally, we define an isotropic objective as the squared Wasserstein distance (E Iso (T ) = 𝒲 (51) One can notice that E Iso is included in 0, 𝜋 2 12 . To have an anisotropic objective included in the same interval we define it as:

E anIso (T ) = 𝜋 2 12 -E Iso (T ) (52)
The gradient of the anisotropic objective is simply the opposite of the isotropy gradient (51).

Then we regroup these two objectives inside the orientation objective defined as: As for the CVT energy, the intersection of T with zone boundaries creates new vertices. Again the gradient is transferred to the initial vertices of T using the chain rule (10).

Finally, the last weight term 𝛾 Obj allows us to control the influence of E Ali and E Ori . We use the following normalized weight: (54)

Zoning objective

The zoning objective is only used during the combinatorial optimization. It counts the number of intersections between T and material zone boundaries:

E Zon = ∑ (𝑖,

  [LXM * 19] explore space-filling curves to obtain close to isotropic mechanical behaviors in the produced part -while typical raster fill patterns (zig-zags) lead to anisotropic mechanical responses. Papacharalampopoulos et al. [PBS18] explore the use of Hilbert curves for continuous deposition, noting that uninterrupted extrusion reduces print time. Continuity is also an essential factor in final quality when extruding pastes or ceramics [HHLT19]. Giannatsis et al. [GVCD15] construct space-filling curves following a material gradient. Chen et al. [CSG * 17] vary the local density of the trajectories following an input image. Zhao et al. [ZGH * 16]

Figure 2 :

 2 Figure2: General overview of our approach. We start by covering the input polygon P with a graph and finding an initial cycle. Then, we proceed to optimize the topology of the graph. Finally, we optimize its geometry. Each frame visually represents a step of the method, and is related to a different section.

Figure 3 :

 3 Figure 3: Initial seeds (in green) are created by intersecting a triangular grid (side length 𝑑) with P.

Figure 4 :

 4 Figure 4: Obtaining aligned edges with the needles. Top: vector field: the green, yellow, and red zones have alignment vectors specified. Middle: standard point Voronoi tessellation. Bottom: Voronoi tessellation obtained using needles. Many of the diagram edges now follow the input vector field.

  submitted to COMPUTER GRAPHICS Forum (4/2022).

Figure 5 :

 5 Figure 5: Closeup on two cells 𝐴 and 𝐵 of a three connected, bridge-less planar graph. Left: a perfect matching of edges (thick black lines). Right: inverting the perfect matching creates cycles covering all graph nodes. Here two cycles (in green and orange color) are produced from which we see only a few edges.

Figure 6 :

 6 Figure6: Top left: the insertion of nodes 𝑛 1 , 𝑛 2 , 𝑛 3 , and 𝑛 4 (in red color) produces new triangular faces. We seek to reconnect the new nodes with the two cycles passing around (in green and yellow color) while covering all nodes. Top right: we remove edges 𝑒 6 , 𝑒 10 and 𝑒 11 and reconnect their endpoints with 𝑛 1 , 𝑛 3 and 𝑛 4 , respectively. However, we need to find an alternative since the resulting cycles do not cover the node 𝑛 2 . Bottom: here we remove edges 𝑒 5 and 𝑒 6 and reconnect their endpoints with nodes 𝑛 1 and 𝑛 2 respectively. This leads to a satisfactory result covering all nodes.

Figure 7 :

 7 Figure 7: Connecting cycles through split-and-reconnect. Here we have two different cycles (in green and yellow), and we identify an opportunity to split-and-reconnect using edges 𝑒 1 , 𝑒 2 , 𝑒 3 and 𝑒 4 (inside red quad). Left: initial configuration, the edge 𝑒 1 belongs to the yellow cycle, and the edge 𝑒 2 belongs to the green cycle. Right: the two cycles are joined by removing 𝑒 1 , 𝑒 2 and adding 𝑒 3 , 𝑒 4 .

  submitted to COMPUTER GRAPHICS Forum (4/2022).

Figure 8 :

 8 Figure 8: Connecting cycles through edge doubling. Left: an unused edge 𝑒 1 (in red quad) between two different cycles (in green and yellow). Right: the two cycles are joined together after doubling the edge 𝑒 1 and its endpoints, adding an extra edge 𝑒 1𝑏𝑖𝑠 .

Figure 9 :

 9 Figure9: Illustration of the four local operations evolving an H-cycle while preserving its Hamiltonian property. Thick green edges belong to the cycle; solid black edges do not belong to the cycle. The state of other edges (gray dashed) has no influence. Operation 1: two triangles given by the edges (𝑒 1 , 𝑒 2 , 𝑒 3 ) and (𝑒 4 , 𝑒 5 , 𝑒 6 ) share one node. We change the cycle edges from (𝑒 2 , 𝑒 3 , 𝑒 6 ) to (𝑒 1 , 𝑒 4 , 𝑒 5 ). Operation 2: two triangles given by the edges (𝑒 1 , 𝑒 3 , 𝑒 4 ) and (𝑒 2 , 𝑒 3 , 𝑒 5 ) share one edge 𝑒 3 . We change the cycle edges from (𝑒 4 , 𝑒 3 , 𝑒 2 ) to (𝑒 1 , 𝑒 3 , 𝑒 5 ). Operation 3: two rhombuses enclosed by edges (𝑒 1 , 𝑒 2 , 𝑒 3 , 𝑒 4 ) and (𝑒 3 , 𝑒 5 , 𝑒 6 , 𝑒 7 ) share a common edge 𝑒 3 . We change the cycle edges from (𝑒 2 , 𝑒 3 , 𝑒 4 , 𝑒 6 ) to (𝑒 1 , 𝑒 3 , 𝑒 5 , 𝑒 7 ). Operation 4: two rhombuses enclosed by edges (𝑒 1 , 𝑒 2 , 𝑒 3 , 𝑒 4 ) and (𝑒 5 , 𝑒 6 , 𝑒 7 , 𝑒 8 ) share a common node and there exists a path in the green cycle (indicated by the green dashed curve) that connects two nodes and that does not traverse any other node in the illustration. We change the cycle edges from (𝑒 1 , 𝑒 4 , 𝑒 5 , 𝑒 7 ) to (𝑒 2 , 𝑒 3 , 𝑒 6 , 𝑒 8 ).

Figure 10 :

 10 Figure 10: Result of our geometry optimization over the brain model with an isotropy objective.

  submitted to COMPUTER GRAPHICS Forum (4/2022).

Figure 11 :

 11 Figure 11: Two-stage optimization of shape depicting a brain. The trajectory of the left part is under an alignment objective, while the trajectory in the right part follows an isotropy objective. The plots depict the evolution of the value of each objective function during the optimization steps. Note that some values are multiplied by a constant factor in order to improve the visualization.

Figure 12 :

 12 Figure 12: Radial vector field in a disc. Left: Our result. Right: Without any global optimization of orientation (initial graph construction directly followed by geometric optimizer without E Ali , E Ori ). The result no longer closely follows the desired orientation.

Figure 13 :Figure 14 :

 1314 Figure13: Correlation between 𝑑 (controlling the seed densitysee Figure3) and the interspace 𝜏 around the optimized cycle. The thick curve is the mean observed over 30 models using different objectives. The top and bottom curves are the reported min/max values of the mean across models.

  submitted to COMPUTER GRAPHICS Forum (4/2022).

𝛾 Lap = 0 𝛾 Lap ≠ 0 Figure 15 : 2 𝑑

 00152 Figure 15: Optimization results showing the effect of using or not the smoothing term 𝛾 Lap .

  submitted to COMPUTER GRAPHICS Forum (4/2022).

Figure 17 :

 17 Figure 17: Left: input specifying zones and objective selection in each. Parallel vec field means that we align the path with a vector field parallel to the zone border and ortho vec field means that we use alignment with a vector field orthogonal to the border. Middle: cycle optimized for a width of 1.7𝑚𝑚. Right: cycle optimized for a width of 0.8𝑚𝑚.

Figure 18 :

 18 Figure 18: Left: Single material leaf printed with a constant width smaller than nominal to appreciate the cycle folds. Note how the cycle goes around holes with a single folded trajectory in the center part. Middle: Multi-material version using a narrower spacing, with the inner skeleton in a different material. Right: Closeup. Note how the cycle locally has a single trajectory around the holes. Such cycles cannot be obtained by growing from a small circle to invade the surface: this always gives an even number of paths in between outer contours.

  submitted to COMPUTER GRAPHICS Forum (4/2022).

Figure 19 :

 19 Figure 19: Multi-material printing of the Batman logo following the parallel vector field.

Figure 20 :

 20 Figure 20: Zebra printed with black and white filaments.

Figure 24

 24 Figure 24 compares extrusion quality between a zigzag infill and our technique in a single layer, printed with flexible filament. The travel moves required by the standard zigzag produce artifacts, as the print head leaves marks and material oozes along with some travel moves. In contrast, our cyclic infill produces a similarly oriented deposition without travel moves. Other examples of near zigzag anisotropic fills are visible in Figure 1, (b).Figure 22 also compares our infill to a standard zigzag on full 3D prints. Figure 25 provides a visual comparison to [ZGH * 16] and [PS06], two techniques producing cyclic trajectories.

  Figure 24 compares extrusion quality between a zigzag infill and our technique in a single layer, printed with flexible filament. The travel moves required by the standard zigzag produce artifacts, as the print head leaves marks and material oozes along with some travel moves. In contrast, our cyclic infill produces a similarly oriented deposition without travel moves. Other examples of near zigzag anisotropic fills are visible in Figure 1, (b).Figure 22 also compares our infill to a standard zigzag on full 3D prints. Figure 25 provides a visual comparison to [ZGH * 16] and [PS06], two techniques producing cyclic trajectories.

  submitted to COMPUTER GRAPHICS Forum (4/2022).

Figure 21 :

 21 Figure 21: (a): plate printed using the per-segment area. (b): plate printed with a constant flow all along the cycle. Both use the same total volume of filament, are printed at the same bed location. (c) top: closeup on the leaf model shown in Figure 18. (c) bottom: closeup on the same leaf printed without the per-segment area.

Figure 23 :

 23 Figure 23: Multi-material cube. Top: Views of the print obtained with a mixing extruder, the transition occurs within the part. Note the isotropic infill pattern on the top. Bottom-left: The white zone has to be traversed when going from blue (mixture A) to red (mixture B), and is where mixture transition occurs. The user can freely place it within each layer. Bottom-right: same model using a zigzag infill, printed on a two-extruders Ultimaker 3. Note the prime tower, requiring extra movements at each layer.

Figure 24 :

 24 Figure 24: Comparing a standard zigzag strategy and our anisotropic infill. Left: While the zigzag is perfectly regular, it does not allow for continuous extrusion. Travel moves create marks and oozing, that damage the final part. Right: None of these issues are present on our fully continuous 3D prints.

Figure 25 :

 25 Figure 25: Left: Result from our technique, with different objectives: vector field orthogonal to border (green), isotropic (orange) and concentric vector field (gray). Note the single crossing near the arms and the feet. The gray part is the one optimized to match the orientation of the two other methods. Middle and right: Results from respectively [ZGH * 16] and [PS06] (images from papers).

Figure 26 :

 26 Figure 26: Comparison between a growth process (a,c) and our approach (b,d). Please refer to the text for details.

Figure 27 :

 27 Figure 27: A segment cell in blue and two point cells in red.

•𝑓Figure 28 :

 28 Figure 28: The segment cell boundary is represented by two functions 𝑓 , 𝑔.

) 2 𝒲𝑊 , 𝑦 + 1 2 ℋ 1 0∫Figure 29 : 1 0∫

 221291 Figure 29: The vector field's integral under a segment is obtained by decomposing it in smaller segments contained in single quads of the grid.

1 0( 2 𝑑𝑡

 12 ) ∈T 𝐸 𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 𝛿 ∠(𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 ) (42)Where ∠𝑣 is the angle of 𝑣 modulo 𝜋 taken inside the interval [0, 𝜋[, and 𝐿 is the total length of T . We know from[START_REF] Rabin | Transportation distances on the circle[END_REF] that the Wasserstein distance on a circle can be expressed as:𝒲 2 2 (U, D T ) = inf 𝛼 ∫ 𝐹 𝛼 ) (-1) (𝑡) -𝐺 (-1) T (𝑡)(43)submitted to COMPUTER GRAPHICS Forum (4/2022).where 𝐹 and 𝐺 T are, respectively, the cumulative distribution of U and D T . That is to say:𝐹 (𝜃) = 𝜃/𝜋, 𝐺 T (𝜃) = 1 𝐿 ∑ (𝑖 𝑗) ∈ T 𝐸𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 𝟙 {∠(𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 )<𝜃} (44)

  ) ∈ T 𝐸 ∠ 𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 𝟙 𝑡 ∈ 𝑇 0 𝑖 𝑗 ,𝑇 1 𝐸 ∠ (𝑣 𝑣 𝑣 𝑙 -𝑣 𝑣 𝑣 𝑘 )<∠(𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 ) 𝑣 𝑣 𝑣 𝑙 -𝑣 𝑣 𝑣 𝑘 , 𝑇 1 𝑖 𝑗 = 1 𝐿 ∑ (𝑘𝑙) ∈T 𝐸 ∠ (𝑣 𝑣 𝑣 𝑙 -𝑣 𝑣 𝑣 𝑘 ) ∠ (𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 )𝑣 𝑣 𝑣 𝑙 -𝑣 𝑣 𝑣 𝑘 (46)

  + 𝜋𝛼 0 -∠ 𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 2 𝑑𝑡 (48)

2 2 ( 2 - 2 - 2 -

 2222 U, D T )). Computing the integral in (48), we get:E Iso (T ) = 1 3𝜋 ∑ (𝑖 𝑗) ∈ T 𝐸 Θ 1 𝑖 𝑗 -∠ 𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖 3 -Θ 0 𝑖 𝑗 -∠ 𝑣 𝑣 𝑣 𝑗 -𝑣 𝑣 𝑣 𝑖3(50) As for the alignment objective, we only compute an approximation of the gradient by ignoring the variations of Θ 𝑘 𝑖 𝑗 . Our approximation of the gradient with respect to 𝑣 𝑣 𝑣 𝑖 is as follows:G (𝑖) Iso = 𝑅 𝜋/2 𝑣 𝑣 𝑣 𝑖 + -𝑣 𝑣 𝑣 𝑖 𝜋 𝑣 𝑣 𝑣 𝑖 + -𝑣 𝑣 𝑣 𝑖 2 𝜋Θ 1 𝑖𝑖 + -∠ 𝑣 𝑣 𝑣 𝑖 + -𝑣 𝑣 𝑣 𝑖 𝜋Θ 0 𝑖𝑖 + -∠ 𝑣 𝑣 𝑣 𝑖 + -𝑣 𝑣 𝑣 𝑖 𝑅 𝜋/2 𝑣 𝑣 𝑣 𝑖 -𝑣 𝑣 𝑣 𝑖 - 𝜋 𝑣 𝑣 𝑣 𝑖 -𝑣 𝑣 𝑣 𝑖 - 2 𝜋Θ 1 𝑖 -𝑖 -∠ 𝑣 𝑣 𝑣 𝑖 -𝑣 𝑣 𝑣 𝑖 - 𝜋Θ 0 𝑖 -𝑖 -∠ 𝑣 𝑣 𝑣 𝑖 -𝑣 𝑣 𝑣 𝑖 -2

EFigure 30 :

 30 Figure 30: Left: A path with 4 segments. Right: Transport map of the path orientation distribution to the uniform distribution. Sticks are Dirac located at the orientation of the segment having the same color. They are mapped to domains of [0, 𝜋[ represented by portions of the ring having an area proportional to the segment length.

  ) 𝑑𝑥 + ∑ 𝑖 area (I 𝑖 ) + ∑ 𝑖 area (A 𝑖 ) -1

  𝑗) ∈T 𝐸 ∑ 𝑘≠𝑙 # [𝑣 𝑣 𝑣 𝑖 , 𝑣 𝑣 𝑣 𝑗 ] ∩ Z 𝑘 ∩ Z 𝑙 (55) Where #𝑆 is the cardinal of a set 𝑆. The weight 𝛾 Comb attempts to prioritize the zoning objective E Zon . To do so, we set 𝛾 Comb such that 𝛾 Comb (E Ali + E Ori ) is included in [0, 1]. The squared angle difference in the alignment objective is bounded by 𝜋 2 /4, and the squared Wasserstein distance is bounded by 𝜋 2 /12. We then set 𝛾 Comb to:

	𝛾 Comb =	𝜋 2 4	∫	-1
				(56)

P 𝑉 (𝑥) 𝑑𝑥 + ∑ 𝑖 area (I 𝑖 ) + ∑ 𝑖 area (A 𝑖 ) submitted to COMPUTER GRAPHICS Forum (4/2022).

submitted to COMPUTER GRAPHICS Forum (4/2022).

Acknowledgements

This research was initiated at the 18th Bellairs Workshop on Computational Geometry, February 16-22, 2019, co-organized by S. Lazard and S. Whitesides. We thank the other workshop participants for stimulating conversations and Pierre-Alexandre Hugron for his help with 3D printing. This work was partly supported by Lorraine Université d'Excellence, ANR-15-IDEX-04-LUE, ANR-17-CE10-0002, and NSERC Discovery Grant.

Appendix A: Appendix

CVT objective

We perform a CVT optimization of a diagram of segments. Let 𝑥 ∈ R 2 be a point and 𝑆 ⊂ R 2 be a bounded set. The infimum Euclidean distance between 𝑥 and 𝑆 is:

In order to cover the entire polygon P without keeping empty areas, one could consider minimizing the following objective:

Minimizing the above objective could result in new crossings between the cycle and material zone boundaries, undermining the zoning objective. To avoid this issue, we minimize the 𝐶𝑉𝑇 objective independently for each material zone:

Computing the intersection of T with Z 𝑖 leads to new vertices at the intersections between segments of T and zone boundaries. Consider an intersecting vertex 𝑢 laying inside a segment ]𝑣 𝑣 𝑣 𝑖 , 𝑣 𝑣 𝑣 𝑗 [. Then there exists a real value 𝑡 ∈]0, 1[ such that: