Robotic manipulation of cloth: mechanical modeling and perception

ROBOTIC MANIPULATION OF CLOTH Mechanical Modeling and Perception franco coltraro Doctoral Programme in Applied Mathematics Facultat de Matemàtiques i Estadística Universitat Politècnica de Catalunya supervisors: Jaume Amorós Maria Alberich-Carramiñana February 2023 Franco Coltraro: Robotic Manipulation of Cloth, Mechanical Modeling and Perception, © February 2023 A thesis submitted to the Universitat Politècnica de Catalunya for the degree of Doctor of Philosophy doctoral programme: Applied Mathematics location: Institut de Robòtica i Informàtica Industrial, CSIC-UPC Barcelona, Spain funding: This project was developed in the context of the project CLOTHILDE (”CLOTH manIpulation Learning from DEmonstrations”) which has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 741930) and was also supported by the Spanish State Research Agency through the María de Maeztu Seal of Excellence to IRI (MDM-2016-0656). supervisors: Jaume Amorós Maria Alberich-Carramiñana A mi Nonna. ABSTRACT In this work, we study various mathematical problems arising from the robotic manipulation of cloth. First, we develop a locking-free continuous model for the physical simulation of inextensible textiles. We present a novel finite element discretization of our inextensibility constraints which results in a unified treatment of triangle and quadrilateral meshings of the cloth. Next, we explain how to incorporate contacts, self-collisions and friction into the equations of motion, so that frictional forces and inextensibility and collision constraints may be integrated implicitly and without any decoupling. We develop an efficient active-set solver tailored to our non-linear problem which takes into account past active constraints to accelerate the resolution of unresolved contacts and moreover can be initialized from any nonnecessarily feasible point. Then, we embark ourselves on the empirical validation of the developed model. We record in a laboratory setting –with depth cameras and motion capture systems– the motions of seven types of textiles (including e.g. cotton, denim and polyester) of various sizes and at different speeds and end up with more than 80 recordings. The scenarios considered are all dynamic and involve rapid shaking and twisting of the textiles, collisions with frictional objects and even strong hits with a long stick. We then, compare the recorded textiles with the simulations given by our inextensible model and find that on average the mean error is of the order of 1 cm even for the largest sizes (DIN A2) and the most challenging scenarios. Furthermore, we also tackle other relevant problems to robotic cloth manipulation such as cloth perception and classification of its states. We present a reconstruction algorithm based on Morse theory that proceeds directly from a point-cloud to obtain a cellular decomposition of a surface with or without boundary: the results are a piecewise parametrization of the cloth surface as a union of Morse cells. From the cellular decomposition, the topology of the surface can be then deduced immediately. Finally, we study the configuration space of a piece of cloth: since the original state of a piece of cloth is flat, the set of possible states under the inextensible assumption is the set of developable surfaces isometric to a fixed one. We prove that a generic simple, closed, piecewise regular curve in space can be the boundary of only finitely many developable surfaces with nonvanishing mean curvature. Inspired by this result we introduce the dGLI cloth coordinates, a low-dimensional representation of the state of a piece of cloth based on a directional derivative of the Gauss Linking Integral. These coordinates –computed from the position of the cloth’s boundary– allow us to distinguish key qualitative changes in folding sequences. v RESUMEN En este trabajo estudiamos varios problemas matemáticos relacionados con la manipulación robótica de textiles. En primer lugar, desarrollamos un modelo continuo libre de locking para la simulación física de textiles inextensibles. Presentamos una novedosa discretización usando elementos finitos de nuestras restricciones de inextensibilidad resultando en un tratamiento unificado de mallados triangulares y cuadrangulares de la tela. A continuación, explicamos cómo incorporar contactos, autocolisiones y fricción en las ecuaciones de movimiento, de modo que las fuerzas de fricción y las restricciones de inextensibilidad y colisiones puedan integrarse implícitamente y sin ningún desacoplamiento. Desarrollamos un solver de tipo conjunto-activo adaptado a nuestro problema no lineal que tiene en cuenta las restricciones activas pasadas para acelerar la resolución de nuevos contactos y, además, puede inicializarse desde cualquier punto no necesariamente factible. Posteriormente, nos embarcamos en la validación empírica del modelo desarrollado. Grabamos en un entorno de laboratorio -con cámaras de profundidad y sistemas de captura de movimientolos movimientos de siete tipos de textiles (entre los que se incluyen, por ejemplo, algodón, tela vaquera y poliéster) de varios tamaños y a diferentes velocidades, terminando con más de 80 grabaciones. Los escenarios considerados son todos dinámicos e implican sacudidas y torsiones rápidas de los textiles, colisiones con obstáculos e incluso golpes con una varilla cilíndrica. Finalmente, comparamos las grabaciones con las simulaciones dadas por nuestro modelo inextensible, y encontramos que, de media, el error es del orden de 1 cm incluso para las telas más grandes (DIN A2) y los escenarios más complicados. Además, también abordamos otros problemas relevantes para la manipulación robótica de telas, como son la percepción y la clasificación de sus estados. Presentamos un algoritmo de reconstrucción basado en la teoría de Morse que partiendo directamente de una nube de puntos obtiene una descomposición celular de una superficie con o sin borde: los resultados son una parametrización a trozos de la superficie de la tela como una unión de celdas de Morse. A partir de la descomposición celular puede deducirse inmediatamente la topología de la superficie. Por último, estudiamos el espacio de configuración de un trozo de tela: dado que el estado original de la tela es plano, el conjunto de estados posibles bajo la hipótesis de inextensibilidad es el conjunto de superficies desarrollables isométricas a una fija. Demostramos que una curva genérica simple, cerrada y regular a trozos en el espacio puede ser el borde de un número finito de superficies desarrollables con curvatura media no nula. Inspirándonos en este vii resultado, introducimos las coordenadas dGLI, una representación de dimensión baja del estado de un pedazo de tela basada en una derivada direccional de la integral de enlazamiento de Gauss. Estas coordenadas -calculadas a partir de la posición del borde de la telapermiten distinguir cambios cualitativos clave en distintas secuencias de plegado. viii AGRADECIMIENTOS A mis tutores Jaume y Maria por hacer de guías en lo académico y muchas veces en lo personal en este largo camino. A Amanda por su apoyo emocional. A mis padres y hermano por su amor incondicional. A mis amigos de Barcelona y Madrid por sus palabras de ánimo cuando más las necesitaba. A mis estudiantes de TFG y TFM por todo lo que aprendimos juntos. Al personal del IRI y en especial al Grupo de Percepción y Manipulación por ayudarme, especialmente en mis peleas con el WAM. A todos ellos, Gracias. ix xvi list of figures Figure 3.4 Simulated final frame of cloth’s collision with a collection of needle-like obstacles seen from 3 different angles. The cusps must be taken into account separately from the rest of the surface and are treated with the same algorithm we treat self-collisions. 55 Figure 3.5 Trajectory of the controlled nodes for the dynamic folding of shorts. 56 Figure 3.6 Simulated sequence of the dynamical folding of a pair of shorts. The first part of the motion (frames two and three) is performed fast enough so that the shorts lay partially flat on top of the table after a lowering phase (frame four). The final fold is completed by dropping the top two corners on top of the leg loops (frames five and six). 57 Figure 4.1 In the picture we can see all the fabrics (size A3) used in the experiments. From left to right we have: paper, polyester, light cotton, felt, wool, denim and stiff cotton. 60 Figure 4.2 Experimental setup for the recording of the motion of the real textiles. On the left, the depth camera. On the right, the robotic arm with an affixed hanger. 61 Figure 4.3 Point-cloud obtained using a depth camera (bottom) and RGB image (top). Note that the camera only gets the depth of the objects visible to it, all the rest (e.g. the robot behind the cloth) is occluded. 62 Figure 4.4 Quadrilateral meshing (left) of one of the frames of the filtered and de-noised point-cloud (right). Each quadrilateral is divided into triangles for plotting purposes. 63 Figure 4.5 Comparison at four time instants of the recorded fast motion of paper (right) versus its simulation with the inextensible model (left) with δ=0.37 and α=2.49 . The mean absolute error is 0.41 cm. 64 Figure 4.6 Mean of the absolute error (4.3) for the fast motion of paper using different values of the parameters of the model (damping α and virtual mass δ ). In red we highlight the error with the optimal value of the parameters. 65 list of figures xvii Figure 4.7 Comparison at 4time instants of the fast motion of light cotton (right) versus its simulation with the inextensible model (left) with δ=0.52 and α=2.69 . The mean absolute error is 0.31 cm. 66 Figure 4.8 Comparison of the fast motion of light cotton between the inextensible and the other 3models’ errors (bottom curves) and dispersions (top curves), with respect to the recorded motion (using a 9×9 meshing). The mean absolute error for the inextensible model is 0.31cm. 67 Figure 4.9 Comparison at 4time instants of the slow motion of wool (right) versus its simulation with the inextensible model (left) with δ=0.47 and α=1.19 . The mean absolute error is 0.37 cm. 68 Figure 4.10 Absolute error (bottom curves) and dispersion (top curves) of the position of the simulated textile vs. the recorded motion (with the depthcamera and WAM robot) for the slow movement (top) and the fast one (bottom) using the inextensible model. 69 Figure 4.11 Reflective markers attached to the denim sample (encircled in red). The markers are very small, with a diameter of 3 mm and a weight of 0.013 g. We use 20 reflective markers when the size of the textile is A2and 12 when it is A3.70 Figure 4.12 Setup used to record the motion of the textiles: 5cameras surround the scene so that every marker (highlighted in red in the photo) is visible to at least 2cameras at the same time. This ensures that the system can be certain of the 3D position of the marker. 71 Figure 4.13 Shaking motion sequence (left to right): the cloth is shaken back and forwards. 72 Figure 4.14 Twisting motion sequence (left to right): the cloth is rotated with respect to the z -axis back and forth several times. 72 Figure 4.15 Surface plot of the error function ¯ e(δ,α) for the fast shaking motion of A2wool (left) and a close-up near the detected minimum (right). Notice the presence of noise in the close-up. 74 xviii list of figures Figure 4.16 Three frames comparing the recorded fast twisting of A2polyester (left) with its inextensible simulation (right). The error at the three depicted frames from left to right is 1.82,1.73 and 1.36 cm respectively; being the average error of the whole simulation 1.15 cm. 75 Figure 4.17 Three frames comparing the recorded fast shaking of A2denim (left) with its inextensible simulation (right). The error at the three depicted frames from left to right is 0.63,0.92 and 1.001 cm respectively; being the average error of the whole simulation 0.84 cm. 75 Figure 4.18 Putting a tablecloth motion sequence (right to left): the cloth starts suspended and is afterwards laid dynamically (only partially) onto the table. 78 Figure 4.19 Three frames comparing the recorded tablecloth low friction scenario of A2wool (left) with its inextensible simulation (right). The error at the three depicted frames from right to left is 0.80,1.17 and 0.76 cm respectively; being the average error of the whole simulation 0.58 cm. 79 Figure 4.20 Sensitivity analysis for high friction case, i.e. we vary the value of µ and compute the absolute error (4.3) for the four A2fabrics. In red we encircle the error found with the optimal parameter of µ.79 Figure 4.21 Long-stick hits sequence (left to right): the cloth is held by its two upper corners and then is hit repeatedly with a long stick. The hits are aimed at different locations with varied intensities. 80 Figure 4.22 Long-stick (with a length of 75 cm and a diameter of 1.5cm) used to hit the textiles. Two markers are put at both ends of the stick to record its trajectory. 81 Figure 4.23 Four frames comparing the recorded hitting of A2stiff-cotton (left) with its inextensible simulation (right); being its average error 1.07 cm. On the right, we show a full plot (vertically) of the absolute error and with yellow lines, we highlight the moments in which the stick is in contact with the cloth. 83 Figure 5.1 Critical points of the Morse-Smale function f(x,y,z) = zon an example surface. 88 list of figures xix Figure 5.2 Three types of critical points (from top to bottom: minima, saddles and maxima) and their Morse data for surfaces without boundary. 90 Figure 5.3 Four types of critical points located at the boundary and their tangential and normal Morse data. 92 Figure 5.4 Change in level set when crossing a critical value in a surface (left) and point cloud (right): note the change in neighbors among the 4marked points after the flow. 95 Figure 5.5 Typical problems associated to k-nearest neighbors (left) and Voronoi neighbors (right). On the left, the majority of the closest points to v are clustered at one side of it. On the right, vertices that are too far apart from each other have a neighboring cell. 96 Figure 5.6 In principle, a boundary point can be identified easily because after projecting it and its neighbors on the tangent plane, they cluster in a semi-space of R2 (left panel). Nevertheless, this is not always the case for every boundary point (middle panel). To overcome this, we declare a point as a boundary point when none of the plane projections of its neighbors enclose the point (right panel). 98 Figure 5.7 Intersection of the oriented graph Gdown with a plane. Changes in the number of connected components of the reconstructed curves Γ(c) = Gdown ∩Hc≈ ∪ S1tell us that we have crossed critical points. 100 Figure 5.8 Apparition of a maximum when following the downwards flow: a new connected component appears. 101 Figure 5.9 Disappearance of a minimum when following the downwards flow: an existing connected component disappears. 101 Figure 5.10 Apparition of a saddle point when following the downwards flow: a connected component appears or disappears. 101 Figure 5.11 Different generic and local level set transformations for surfaces without boundary. The reconstructed curves Γ(c) are homeomorphic to a disjoint union of S1.102 xx list of figures Figure 5.12 Different local level set transformations for all possible cases involving boundary points. The reconstructed curves Γ(c) are homeomorphic to a union of S1 and closed intervals [0, 1] . The previous cases (interior minima, maxima and saddles) are not shown again but can also occur. Saddle points are also drawn, in blue. 103 Figure 5.13 A sampled dumbbell: the black line is the direction of the height function; local maxima, resp. minima, are painted red, resp. black; saddle points are painted blue; their 1–cells are outlined in blue. There are two 2 -cells: one in dark red (left) and one in light blue (right). 104 Figure 5.14 Parametrization of a 2 -cell (left) using a rectangle D with isometric sides to the boundary of the 2-cell (right). 108 Figure 5.15 A sampled knotted torus: the black line is the direction of the height function; local maxima, resp. minima, are painted red, resp. black; saddle points are painted blue; their 1 -cells are outlined in blue. On the bottom, we plot the level set curves highlighting when critical points appear. 110 Figure 5.16 A sampled vest: the black line is the direction of the height function; maxima are painted in red, minima in black; 1 -cells corresponding to boundary minima are outlined in blue and the boundary curves in black. The two purple points where the 1 -cells meet each other or the boundary curves are added to the decomposition. The numbers correspond to the different formal 1 -cells that, when identified (e.g. 7with 7’), reconstruct the entire surface from 2pieces homeomorphic to disks. 112 Figure 5.17 Parametrization by a rectangle of the rightmost 2 -cell of Figure 5.16. The bounding 1 -cells are mapped isometrically to a flat rectangle and then interior points are obtained using the neighboring relationships of the cloud. The 2 - cell consisting of 27 000 points (shown in red) is down-sampled interpolating linearly to a cloud of 900 points (shown in blue). 113 list of figures xxi Figure 5.18 A3D-scan of real pants: the black line is the direction of the height function; maxima are painted red, minima in black, saddle points are painted blue; 1 -cells are outlined in blue and the boundary curves in black. The purple point where the 1 -cells meet is added to the decomposition. On the bottom, we plot the level-set curves obtained by intersecting the surface with planes perpendicular to the height function. 115 Figure 6.1 A smooth simple closed curve (in black) which is the boundary of two developable surfaces (indicated in red and blue respectively) 123 Figure 7.1 Folding sequence of a quadrangular cloth with its associated dGLI cloth coordinates, represented as upper triangular matrices. Each matrix element mij is a geometrical value corresponding to the dGLI between the segments i and j highlighted in red in the corresponding folded state. Notice how some values of the matrix change sign when corners are folded or cross each other. 126 Figure 7.2 The subset of chosen segments is marked in red. 132 Figure 7.3 Study of the index during 3folding sequences. In the left column, we show a representation of the cloth frames, and in the right column the confusion matrix of all of them. In red we highlight the clear class changes that can be identified. 134 Figure 7.4 Confusion matrix that computes all the distances between the states shown in the top table. 135 Figure 7.5 Synthetic representatives chosen for each class. When only one is chosen, it is the closest to the centroid of the class. When a class has more sparsity, additional representatives are chosen to represent the subgroups in the class. 139 Figure 7.6 Results of the real image classification using the data-base presented in Figure 7.4as reference. The first column shows the ground truth class of the images, and at the bottom of every image the classified class. 140 Figure A.1 Four frames comparing the recorded hitting of A2denim (left) with its inextensible simulation (right); being its average error 0.98 cm. On the right, we show a full plot (vertically) of the absolute error and with yellow lines we highlight the moments in which the stick is in contact with the cloth. 155 Figure A.2 Four frames comparing the recorded hitting of A2wool (left) with its inextensible simulation (right); being its average error 1.39 cm. On the right, we show a full plot (vertically) of the absolute error and with yellow lines we highlight the moments in which the stick is in contact with the cloth. 156 LIST OF TABLES Table 2.1 Physical parameters of the inextensible model and their meaning. 27 Table 2.2 Distances (cm) between the bottom corner of the cloths (of dimensions 1m×1m). 30 Table 4.1 Density, sizes and examples of all the materials used in all the experiments. 60 Table 4.2 Estimated parameters and mean absolute errors (cm) for the slow and fast oscillatory motions performed by the WAM robot. 64 Table 4.3 Comparison of the 4models. The first column is the mean of the absolute error (4.3) of the simulated cloth with respect to the recorded motion. 67 Table 4.4 Average velocities ( m·s−1 ) for the twisting motion of the A2textiles. We display the speeds for the first repetition (with a hanger) and for the second (with bare hands). 73 Table 4.5 Mean absolute error ¯ e and the mean standard deviation ¯ s with the optimal value of the parameters averaged over: fabric’s material, size (A2or A3), type of movement and speed for the first repetition of the recordings. 74 Table 4.6 Optimal values of the parameters δk,αk for both repetition I (with hanger) and II (with bare hands) obtained by minimizing the function R.76 xxii list of tables xxiii Table 4.7 Optimal values of the friction coefficients along with the mean absolute error and spatial standard deviation for the low and high friction scenario 78 Table 4.8 Mean absolute error and spatial standard deviation with the optimal value of the parameters α∗ and δ∗ . In the two last columns, we display the quotient (4.9), for our active-set collision algorithm and a standard interior-point method. 82 Table 4.9 Comparison of the parameters ( ˆ δ,ˆ α ) obtained with the aerodynamic formulas (4.10) with the optimal ones ( δ∗,α∗ ) along with their respective mean absolute errors. 84 Table 7.1 Comparison between different shape representations* 137 Table A.1 Summary of results of the aerodynamic study for the first repetition (with bare hands) of the experiments. For all the 32 recordings we display the characteristics of the recording (fabric, size, motion and speed), and the optimal value of the fitted parameters along with their associated absolute error and spatial standard deviation. 153 Table A.2 Summary of results of the aerodynamic study for the second repetition (with hanger) of the experiments. For all the 32 recordings we display the characteristics of the recording (fabric, size, motion and speed), and the optimal value of the fitted parameters along with their associated absolute error and spatial standard deviation. 154 1 INTRODUCTION Robotic manipulation of cloth in a domestic environment, i.e. not in its manufacturing plant, is an increasingly relevant problem because of the ubiquity and versatility of textiles in human lives and activities; with promising applications ranging from folding clothes to dressing people with impaired mobility [31,40]. One of the main challenges faced is the high variety of deformation states that textiles can present [100] (see Figure 1.1). In contrast to rigid body manipulations, where the dynamics of the manipulated object are very well understood [110], there is not one single model that can be considered best in terms of describing the dynamics of real textiles [85]. Having a faithful and easy-to-adjust physical model of cloth behavior is useful for planning and control strategies in the task of robotic cloth manipulation [26,71]; as well as for generating the massive data required to train learning algorithms before their deployment and tuning in the real world [25, 58]. Figure 1.1: Various deformation states that a shirt can present: crumpled, extended, hanging, etc. Image adapted from [61]. Therefore, automatic manipulation of textiles is being tackled nowadays by combining analytical modeling with data-driven learning, so as to overcome the shortcomings of either individual approach. Cloth models need to retain the physical properties relevant for dynamic manipulation and, to be time-affordable, dispense with other features such as appealing rendering. The need for a compact, quick-tocompute model for deformable objects –and textiles, in particular– has 1 8 introduction account for friction between cloth and possible obstacles (e.g. a table) and with itself. This is done in a manner that will allow us to consider later all constraints (inextensibility and contacts) and friction forces at the same time without any decoupling. In Section 3.3we explain how to detect and include self-collision constraints under the framework presented in Section 3.2. Particularly important for efficiency and to avoid unwanted oscillations is how to take into account the thickness of cloth. This is explained in Section 3.3.3. Afterwards, in Section 3.4 a novel numerical discretization is presented in order to integrate the extended equations of motion. This can be seen as a natural extension of the fast projection algorithm presented in [44] in order to include inextensibility, contacts and friction in a single pass. This discretization leads naturally to a sequence of quadratic problems with inequality constraints. We explain in detail how to include selfcollision constraints under this scheme. In Section 3.5we enter the most mathematical part of the chapter, where we study how to solve efficiently the sequence of quadratic problems defined before. We present a novel active-set method tailored to our problem. The main advantage of this new algorithm will be its ability to start from any point and not necessarily from a feasible one. Lastly, we discuss how to implement it efficiently with the use of Cholesky factorizations. A detailed procedure is laid out in pseudo-code in Algorithm 2. To close the chapter we present several scenarios that put to a test the developed collision model. They will be qualitative in nature, i.e. we only show that our simulator is capable of dealing with them (more quantitative experiments will be performed in Chapter 4, Section 4.3). We will show that the model of friction is effective in static (cylinder experiment, Section 3.6.1) and dynamic (rotating sphere experiment, Section 3.6.2) settings. Next, we show how we can easily include collisions with sharp needle-like objects in Section 3.6.3. Finally, we simulate complicated folding sequences of cloth with non-trivial topologies (a pair of shorts) in Section 3.6.4. The second and the third experiments are challenging scenarios suggested by [70] as tests that every robust collision model should be able to pass. 1.3.3Chapter 4: experimental validation In Section 4.1we compare the inextensible model with recorded motions of real fabrics under different scenarios. First, we show how with as few as two parameters, we are able to model within an error of less than 0.5cm different types of DIN A3textiles (such as cotton, wool and felt) under fast and slow shaking movements. To perform this real-world validation we use a Barret robotic arm together with a depth camera to record the real motion of garments being shaken at different velocities. Afterwards, in Section 4.1.3we perform a sensitivity analysis in order to understand how stable is our model with 1.3 organization 9 respect to the fitted parameters. Finally, the performance of the inextensible model is compared to other popular cloth models (Section 4.1.4). Figure 1.4: Four frames comparing the recorded hitting of A2polyester (left) with its inextensible simulation (right); its average error being 1.44 cm. On the right, we show a full plot (vertically) of the absolute error and with yellow lines, we highlight the moments in which the stick is in contact with the cloth. In the second round of experiments, we perform a more exhaustive set of recordings than before. We enlarge our database by adding a new twisting movement and a larger cloth size (DIN A2). The fabrics are shaken and twisted by a human at two different speeds and are recorded with a Motion Capture System. We carry out two repetitions of each motion and at the end have 64 different recordings of about 15 10 introduction seconds. As before we estimate the optimal values of the two physical parameters and achieve mean errors of less than 1cm (see Tables 4.5, A.1and A.2), even for the A2textiles and fast motions. Finally, we try to understand how the speed and size of the textiles affect their motion by developing a predictive and a priori formula for the value of the cloth’s physical parameters (see Section 4.2.3). To close up the chapter we perform a final set of experiments where we use again the motion capture system, this time to record the collision of four (size A2) textiles in two different scenarios. In one of them, the fabrics are laid dynamically on top of a table in a putting-a-tablecloth fashion. In the other, they are hit by a long stick four times at various places and with different strengths (see Figure 1.4). For the first scenario, we find the optimal friction parameter for both a high and a low friction case and study the stability of the model with respect to this parameter (see Figure 4.20). For the hitting experiment, we find again the optimal parameters of the model but this time we also put to use the predictive formulas found before in Section 4.2.3by using them to compute an estimate of the physical parameters without performing any optimization (see Table 4.9). To conclude, we check computational times and compare our active-set solver with a standard interior-point method. 1.3.4Chapter 5: surface reconstruction via Morse theory In this chapter, we present an algorithm for the reconstruction of surfaces from point samples. Our method proceeds directly from the point-cloud to obtain a cellular decomposition of the surface derived from a Morse-Smale function (see Section 5.2). After running the algorithm we obtain a piecewise parametrization of the surface as a union of a small number of Morse cells, suitable for tasks such as noiseremoval or meshing, and a cell complex of small rank determining the topology of the manifold (see Section 5.4for all the implementation details). We explain how to apply this algorithm to samples of smooth surfaces with or without boundary, embedded in an ambient space of any dimension. 1.3.5Chapter 6: developable surfaces In this short chapter, we study from a geometric viewpoint developable surfaces, i.e. smooth surfaces with Gaussian curvature 0 . The original state of a piece of cloth is flat, so the set of possible states under the inextensible assumption is the set of developable surfaces isometric to a fixed one. Section 6.1discusses two candidates to the role of generalized coordinates in the space of states of a developable surface under our isometric strain model, and explains their common limitation from the viewpoint of their application. Section 6.2proposes 1.4 videos 11 an alternative approach: to track the motion of the surface by following its boundary. This is not straightforward because the boundary does not determine completely the position of the surface, but as we explain below our Main Theorem 4establishes the feasibility of this approach. 1.3.6Chapter 7: semantic classification of cloth states In this chapter, we define the dGLI Cloth Coordinates, a low-dimensional representation of the configuration state of a rectangular cloth piece that allows us to distinguish topologically different folded states. In Section 7.3we introduce the novel concept of the directional derivative of the GLI (Gauss Linking integral) which allows us to overcome the main limitation of the GLI in planar settings: it becomes zero and therefore uninformative. We derive first an expression for the GLI of two segments, then we prove that we can perturb the segments slightly to obtain information when they are co-planar. Finally in Section 7.4, we apply these new coordinates to a data-base of cloth configurations taken from simulated folding sequences and then we test experimentally the index on real images of folded clothes. 1.4 videos Along this work, several videos will be presented in order to illustrate and show several simulations. All videos (in order or appearance) can be found in the playlist: https://youtube.com/playlist?list=PL1XXvX-KsL9puQ79M03K8A hzCoXxNIVul 1.5 notation Before we proceed, for reference, we present a list of the most important symbols used in this thesis. Moreover, for Part I (the most involved notationally) we will follow the following notational conventions: 1. Matrices, tensors and vector functions will be denoted by bold capital letters. 2. Vectors (except for points in R3 ) will be denoted by bold lowercase letters. 3. The rest (e.g. scalar parameters, points in R3 , scalar functions, etc.) will be denoted by capital and lower-case italics. 12 introduction list of most important symbols SSurface used to model the cloth. ΩeElements (e.g. triangles) of the discretized surface. φe(ξ,η)Local parametrization of the element Ωe. ∂ξφe,∂ηφePartial derivatives of the parametrizations. nNumber of vertices/nodes of the discrete surface. pi(t)Coordinates in R3of the node iof the surface. NiIndicator functions, i.e. Ni(pj) = δij. x(t),y(t),z(t)x-coordinates (resp. yand z) of all the nodes of the surface in Rn. φ(t)Position of all nodes together, i.e. (x(t),y(t),z(t))⊺∈R3n. Ek,Fk,GkCoefficients of the metric of the discrete surface at node k. C(φ(t)) Constraint function, e.g. Ci(φ(t)) = Ek(t)−Ek(0). H(φ(t)) ≥0 Set of collision constraints. λ(t),γ(t)Lagrange multipliers associated to the inextensibility and contact constraints respectively. TE,TF,TGTensors used to compute Ek(t),Fk(t),Gk(t). ˙φ(t),¨φ(t)Velocity and acceleration of all nodes. M,K,DMass, stiffness and damping matrices. fµ(˙φ)Friction force. V(˙φ)Unit relative tangent velocities at the points of contact. dt >0 Time step used to discretize the equations of motion. φm,˙φmPosition and velocity of the nodes at tm=m·dt. φm→dt φm+1Detection of self-collisions between two states. W,OWorking and observation set of the active-set solver. QAQ⊺=LL⊺Cholesky decomposition of the symmetric positive definite matrix A. ϕmNodes of the meshed recording of real textiles at tm. emTime-dependent mean absolute error in cm at tm. smTime-dependent spatial standard deviation in cm at tm. Part I MODELING CLOTH In this first part we deduce the inextensible cloth model: under the assumption that cloth is a surface that moves through space whose metric is preserved, we derive novel inextensibility equations that we later discretize using the Finite Element Method. We test the new model with scenarios involving mesh independency, cloth’s area preservation and complex topology simulations. Then we introduce the problem of collision modeling for inextensible cloth. We develop a new active-set algorithm tailored for resolving contacts (including self-collisions), friction and inextensibility in a single pass without any decoupling of contact and inextensibility constraints. We put to a test the developed collision procedure with scenarios involving static and dynamic friction, sharp objects and complex-topology folding sequences. Finally, we compare the inextensible model developed in the previous two chapters with collected experimental data from real textiles. We study how faithfully our model is capable of reproducing recorded textiles under challenging scenarios involving fast motions, friction and even strong hits with a long stick. 2 INEXTENSIBLE CLOTH MODEL After reviewing the state of the art in Section 2.1, we explain how to impose the inextensibility conditions (2.2) as hard constraints using Lagrange multipliers in Section 2.2. We discretize the cloth using the Finite Element Method (FEM) in Section 2.3and derive an original discretization of the inextensibility equations in Section 2.4. There we discuss practical and computational aspects of evaluating the inextensibility constraints. The system of ordinary differential equations modeling the dynamics of the cloth is presented in Section 2.5along with our physical model of aerodynamics. We discuss how to integrate numerically the equations of motion in Section 2.6. Finally, in Section 2.7we test the performance of the presented inextensible model in several challenging scenarios. 2.1 related work The mechanic behavior of cloth has been extensively studied from the engineering and computer graphics viewpoints. The engineering focus has usually been on cloth manufacture and its mechanical resistance, using elasticity theory to examine local properties of the material. The aim in computer graphics has always been the simulation and representation of cloth motion, using mostly mass-spring models for speed of computation and paying attention to global problems such as nontrivial cloth topology and collisions. In both fields, the textiles are usually modeled as two-dimensional surfaces. Since the pioneering work of [111] dozens of different models have been proposed for cloth simulation. Our approach retakes the original idea of [111]: understanding cloth’s internal dynamics as the preservation of the first fundamental form of the surface. For surveys and research problems about cloth modeling and simulation see [11,22,85]. We now review some of the most popular methods for simulating the internal dynamics of cloth, noting that there is not a clear-cut separation between some of them. We deliberately only reference lines of research that model internal dynamics of cloth in an essentially different physical manner. Approaches that develop novel numerical methods to integrate existing physical models will not be considered. textile engineering: The modeling of fabrics from a Structural Mechanics viewpoint has been a topic of interest in Textile Engineering since the start of the industrial manufacture of cloth. See [50] for a survey, and completion in many respects, of this theory, or [56] for 15 16 inextensible cloth model updates. This modeling has been hampered by the complexity of fabric as a material: it is inhomogeneous, even discontinuous, already at a relatively sizable scale; highly anisotropic and capable of a nonlinear response even to relatively modest strains; prone to buckling under very small compression. Because of these complexities, investigations of the explicit mathematical expression of the stress-strain relationships of fabrics start from the elasticity theory of solids, and typically consider fabric as an elastic thin plate, more rarely as a shell. Ignoring its thickness, the resulting surface has 6strain fields that determine its displacement and deformation, and a 6×6 symmetric stiffness matrix giving the constitutive equations that relate stress and strain ([56] p.10). Particular properties of fabric such as orthotropicity, reduce the number of stiffness parameters to take into account and lead to simplifications of this general elastic model (e.g.[50] and [120]), or to the introduction of nonlinear responses in the model [74]. The development of these models has been based foremost on experimentation, using procedures such as the Kawabata Evaluation System (KES) [65], which accurately measures the tensile, shear, pure bending, compression and friction behavior of a cloth sample under stresses that are typical in the cloth manufacturing process, in order to establish the stiffness parameters of the model. The increase in available computing power has made viable, although not for real-time computing, sophistications of this basic model to take into account the fiber structure of yarns and fabrics (see [23] and [60]). The computational complexity of these models opened the way to the development of simpler, descriptive models that would later become of common use in Graphic Computer Science (as we explain bellow), treating fabrics as a viscoelastic material that can be modeled as a mass-spring system [88]. While such descriptive models have been extensively used in fabric simulation, their tenuous connection to reality has made them of limited use for the industrial handling of cloth. Here the thin-plate with linear elasticity models prevail. elastic models: these models, e.g. [34,99] and [10], derive the internal forces of the garment from elasticity theory. They are the practical realization of the textile engineering models previously discussed, aimed at graphical simulations rather than at the study of physical properties of cloth. They are all continuous in nature and have the advantage of being physically based, stable under different meshing of the garment and convergent when the mesh is refined [22]. Mostly, they use finite elements to discretize the equations of motion [17]. They can be expensive to evaluate (especially when using non-linear elasticity or the co-rotational method, necessary for ensuring rotational invariance [99]) and may exhibit locking of triangular 2.2 inextensibility modeling 17 meshes if elasticity is heavily reduced: see [32] and [63] and more in general [8] and [89] for a more detailed discussion. mass-spring systems: these models, e.g. [93] and [21], derive internal forces from spring-like energies. They are cheap to evaluate and very intuitive, but less physically sound: they are mesh dependent, have a lot of (non-physical) parameters to be tuned (see [81] for an evolutionary algorithm to tune them) and do not show a convergent behavior when the mesh is refined [85]. constrained dynamics: these models derive internal forces from explicit conditions that the cloth must satisfy [11]. Most of the methods use some kind of Lagrange multipliers to impose the constraints. They mostly differ in what the conditions are and the algorithm used to impose the constraints. Our method fits in this category. Sometimes they are used as velocity filters that complement the previous methods. There are mainly two kinds: 1. Continuous: in this case the constraints are continuous in nature, e.g. imposing bounds to the strain tensor. Examples of this approach are [76,84,112] and [118]. To our knowledge, all methods use elasticity theory to some extent, and therefore in the limit, when elasticity is reduced, face the same problems commented previously. Our method lies in this category with the important difference that it does not use elasticity, but differential geometry of the surface to derive the constraints. 2. Discrete: in this case constraints are discrete in nature, e.g. preserving the length of the edges of the meshed garment or the area of the elements. They are derived concretely for the mesh at hand. Examples are [32,44,47,63]. Their main drawback is the possible lack of convergence of the methods (as opposed to the continuous case) and the possible introduction of mesh-related artifacts [76]. Nevertheless, they can be fast and handle better the inelastic scenario than the elasticity-based continuous models. others: There are two different recent trends: modeling woven cloth at the fiber level (e.g. [23] and [6]) as opposed to the macroscopic level (i.e. as a surface) and considering cloth as a non-Newtonian fluid using the Material Point Method (MPM) adapted to co-dimensional elasticity [60]. The first method is computationally very intensive and the second while being more efficient depicts cloth as very elastic. 2.2 inextensibility modeling As explained before, our main assumption will be to consider cloth as acontinuous and inextensible two-dimensional surface. Formulating the model at the continuous level is important because it avoids as much as possible mesh dependencies: i.e. without changing the physical 24 inextensible cloth model Algorithm 1Computation of TE,TF,TG 1:m←#nodes of each element ▷(m = 3(triangles) or m = 4(quads)) 2:for element Ωedo 3:φe:=∑4 i=1pe iNi▷(parametrization of Ωe) 4:I←Index(Ωe)▷(node’s indices of Ωee.g. [4, 48, 50, 13]) 5:for Gauss point xland weight wldo 6:El←φξ(xl)·φξ(xl); 7:Fl←φξ(xl)·φη(xl); 8:Gl←φη(xl)·φη(xl); 9:dAl←q|El·Gl−F2 l| · wl; 10:for k,i,j=1, . . . , mdo 11:tkij F+ = Nk(xl)·∂ξNi(xl)·∂ξNj(xl)·dAl; 12:tkij F+ = 1 2Nk(xl)·∂ξNi(xl)·∂ηNj(xl) + ∂ηNi(xl)·∂ξNj(xl)· dAl; 13:tkij G+ = Nk(xl)·∂ηNi(xl)·∂ηNj(xl)·dAl; 14:end for 15:end for 16:for k,i,j=1, . . . , mdo ▷(global assembly of the tensors) 17:TI(k)I(i)I(j) E+ = 1 mI(k)tkij E; 18:TI(k)I(i)I(j) F+ = 1 mI(k)tkij F; 19:TI(k)I(i)I(j) G+ = 1 mI(k)tkij G; 20:end for 21:end for where we have grouped the previous tensors in just one TE,F,G= [TE;TF;TG]∈Rnc×n×n , and of course C0=TE,F,G⊗P(0) . With this definition, computing the Jacobian matrix ∇C∈Rnc×3n is very easy (because the tensors are symmetric in the last two indexes): ∇C(φ(t)) = 2·[TE,F,G⊗x(t),TE,F,G⊗y(t),TE,F,G⊗z(t)], (2.16) recall that φ(t) = (x(t),y(t),z(t))⊺∈R3n. 2.4.5Practical implementation We now list some considerations to take into account in a practical implementation of our method: - As mentioned before, the constraints (2.7) are added to C to impose that all edges of all boundary curves are preserved. Also, we remark again that the tensor constraints are only applied for interior nodes, except in the presence of corners, where we include their indices in the definition of TFin order to avoid shearing there. - The tensor TE,F,G needs to be computed only once at the beginning of the simulation because it is time independent. All the integrals 2.4 constraint function enforcing inextensibility 25 involved are evaluated exactly using standard Gaussian quadratures (see [17] for details) as shown in Algorithm 1. - In order to be efficient we first calculate the Jacobian (according to Equation (2.16)) and then put C(φ) = 1 2∇C(φ)·φ−C0 , so the matrix P(t)is actually never computed. - The tensor TE,F,G is highly sparse due to the fact that most products of the 3 indicator functions (Equation (2.11)) are zero. It is not difficult to see that it has of the order of ∼150nnon-zero elements. 2.4.6Shearing Energy Although we are assuming that we can model cloth as inextensible (reasonable for denim, stiff cotton, felt, etc.) in all directions, it is known that this is not completely true in other cases (silk, light cotton, wool). Especially relevant for some types of textiles is shearing, i.e. stretching in the diagonal direction. Nevertheless, we will show in Chapter 4that the inextensible assumption is still a very realistic and accurate assumption. For a smooth surface S parametrized by φ , a shearing energy can be modeled (see [111]) as S=ks 2ZS⟨φξ,φη⟩2dA, where ks>0 is a constant that controls the amount of shearing allowed. With the ideas previously presented, we can discretize this energy in a very sound and efficient way. Indeed for a discrete surface Swith nodes φ(t), writing CF(φ(t)) = TF⊗P(t)−C0 F∈Rn(2.17) where TF is defined by Equation (2.14), k∈ {1, . . . , n} and C0 F= TF⊗P(0); we can define the fundamental shearing energy as S(φ(t)) = ks 2CF(φ(t))⊺·M·CF(φ(t)) where M is the mass matrix (see Equation 2.9). It is clear that this shearing energy is invariant under rigid motions (because CF is). Computing the Jacobian ∇CF∈Rn×3n as before (see Equation 2.16), the resulting force derived from this energy is Fs(φ(t)) = −∇S(φ(t)) = −ks∇CF(φ(t))⊺·M·CF(φ(t)) ∈R3n. (2.18) The fundamental shearing energy is the discrete version of a continuous energy. This is important because keeping the value ks>0 fixed, as the mesh is refined (n→+∞) we observe a convergent and stable behavior. 26 inextensible cloth model Remark 2.4.3.For defining the fundamental shearing energy, ideas from [10] were used. Nevertheless, their shearing model is very different from ours: unlike theirs, we derive the discrete shearing energy from a continuous integral equation. From a more practical point of view, they obtain only one condition for every triangle, whereas we have one for each node. The integration of this new energy must be performed implicitly (since it is very stiff), and hence we need to compute the gradient of the force (the Hessian of the energy). As mentioned in [10] the full Hessian leads to an indeterminate matrix (which causes many problems when solving linear systems), and hence it is approximated by: ∇Fs(φ) = −∇2S(φ)≃ −ks∇CF(φ)⊺·M· ∇CF(φ). 2.5 equations of motion of the cloth Finally, we can write down the Lagrangian (kinetic minus potential energies) of our discretized surface: L(φ(t), ˙φ(t)) = ρ 2˙φ(t)⊺·M·˙φ(t)−ρg⊺·M·φ(t) −κ 2φ(t)⊺·K·φ(t)−λ(t)⊺·C(φ(t)), where 1.ρ>0 is the density of the cloth (assuming homogeneous mass), M is the augmented mass matrix and g= (0, . . . , 0|0, . . . , 0|g, . . . , g)⊺ where g=9.8m/s2is gravity, 2. the stiffness matrix (we are using the isometric bending model described in [12]) is K=L⊺ML where L is an approximation of the point-wise Laplacian and κ>0 is a bending constant, 3. and finally λ(t) are the Lagrange multipliers ensuring inextensibility (at the interior nodes and boundaries) and other possible positional constraints (e.g. the corners of the textiles to be manipulated) included. Thus, we get as Euler-Lagrange equations [110] the following ODE system:    ρM¨φ=fρ−κKφ−D˙φ− ∇C(φ)⊺λ C(φ) = 0 (2.19) where fρ=−ρMg is the force of gravity and we have added Rayleigh damping: D=αM+βK, where α and β are positive parameters [127]. Notice how our inextensibility assumption reduces greatly the number of physical parameters of the model (we do not have shearing or 2.5 equations of motion of the cloth 27 stretching parameters and their respective dampings). Nevertheless, cloth dynamics can be very complicated and there is one important factor we are not accounting for in the previous equations: air resistance. Although there exist some simplified models [77], aerodynamic forces on a deformable object (i.e. cloth) submerged in a fluid (i.e. air) are difficult to model (see [41] and [73]), especially near the boundaries of cloth, because turbulences appear and the nonrigid response of cloth to such turbulences is way more unpredictable than that of a rigid, even vibrating, body (see [104]). Remark 2.5.1(Shearing forces).In the case we want to allow some shearing of the cloth, we simply introduce the force Fs(φ(t)) defined in Equation (2.18) into the first equation of the ODE system (2.19). Obviously, in that case, we would have one parameter more ks>0 . For the rest of this thesis unless explicitly stated (namely in Sections 3.6.2,3.6.4and Chapter 7) we assume that no shearing is allowed and set κs=0. parameter meaning ρDensity (inertial mass) δVirtual (gravitational) mass κBending/stiffness αDamping of slow oscillations βDamping of fast oscillations Table 2.1: Physical parameters of the inextensible model and their meaning. 2.5.1Modeling of aerodynamics through virtual mass The equivalence principle states that for any object its inertial mass is equal to its gravitational mass, which means that an object in vacuum of inertial mass m would fall freely under the action of a force of magnitude mg , where g is gravity. Nevertheless, experiments show that in the presence of air, the free-falling velocities depend on the (shape of the) object at hand. In our experiments, we have found that aerodynamic effects can be modeled without explicitly including them by allowing inertial and gravitational masses to be different. Although this is not true in the physical world, the consideration of the two masses as virtual in our model allows us to account for all aerodynamic contributions to cloth dynamics (drag, lift, turbulences at the boundaries) with great accuracy. Hence, we put fδ=−δMg , setting δ as a new parameter to be fitted. In Table 2.1we summarize the meaning of all the parameters of our model. 28 inextensible cloth model 2.6 discretization of the equations of motion As usual, we approximate φ(t) and ˙φ(t) with {φ0,φ1, . . . } and {˙φ0, ˙φ1, . . . } , where φn and ˙φn are the position and velocities of the nodes of the mesh at time tn=n·dt and dt >0 is the size of the time step. Using an implicit Euler scheme to integrate Equations (2.19) we obtain:    φn+1=φn+1 0−M−1∇C(φn+1)⊺λn+1 C(φn+1) = 0, (2.20) where φn+1 0 is the unconstrained step (which depends on φn and ˙φn ). Note that at time tn=n·dt the only unknown in previous equations is φn+1(and λn+1). 2.6.1Fast projection algorithm Finding φn+1 can be done solving the following constrained optimization problem:    minφn+1(φn+1−φn+1 0)⊺·M·(φn+1−φn+1 0) s.t. C(φn+1) = 0, (2.21) because Equations (2.20) are the critical (stationary) points of the optimization problem. This is a quadratic program with quadratic constraints which can be solved with Newton’s method. Nevertheless, in order to make it computationally tractable in front of difficulties such as indefinite system matrices, it is approximated by a sequence of quadratic programs with linear constraints. Write φj+1=φj+∆φj+1 and make the approximation C(φj+1) = C(φj+∆φj+1)≃C(φj) + ∇C(φj)∆φj+1, then the sequence is:    min∆φj+1∆φ⊺ j+1·M·∆φj+1 s.t. C(φj) + ∇C(φj)·∆φj+1=0, the initial point being φ0=φn+1 0 . This is called the fast projection algorithm [44]. Each of these quadratic programs (which has a diagonal matrix system M ) can be reduced to solving a linear system and we iterate until the constraints are satisfied to a given relative tolerance (usually 0.1% ). We have found this algorithm very stable and fast for our purposes. Remark 2.6.1.In the presence of an obstacle (e.g. a table) given by an implicit equation H(φ) = 0 with a well-defined outwards normal ∇H(φ) , we can model collisions by including new constraints in every iteration of the fast projection algorithm described earlier, i.e. we solve 2.7 evaluation and results 29 iteratively the following sequence of quadratic programs with linear equality and inequality constraints:          min∆φj+1 1 2∆φ⊺ j+1·M·∆φj+1 C(φj) + ∇C(φj)∆φj+1=0, H(φj) + ∇H(φj)∆φj+1≥0, where we have made the approximation H(φj+1) = H(φj+∆φj+1)≃H(φj) + ∇H(φj)∆φj+1. A theoretical justification together with an algorithm to solve efficiently the previous quadratic problem will be discussed in detail in Chapter 3. 2.7 evaluation and results In this section, we study experimentally several desirable properties of the inextensible cloth model. The section is organized as follows: first, we perform a simple quasi-static test to show the locking-free nature of our model, along with its independence with respect to different meshed topologies (Section 2.7.1). We use different triangle and quadrilateral meshings to prove our point. Second, we show how to simulate Cusick’s test with our model (Section 2.7.2). With this experiment, we show the stability of our model when the mesh is refined. Lastly, we present a scenario with non-trivial cloth topology, where we simulate the motion of a tank-top and check to what extent our theoretical inextensibility assumptions are being met in practice, computing several area errors (Section 2.7.3). 2.7.1Locking test With this experiment, we intend to show the locking-free nature of our inextensible model (no shearing is allowed κs=0 ), and its stability with respect to mesh topology. We fix three corners of a flat sheet of cloth (with added random noise of standard deviation 3mm) of 1m by 1m, and let the fourth corner fall freely. Figure 2.4: Different meshes used for the locking test. The triangular meshes (left) are irregular whereas the quadrilateral ones (right) are uniform. 30 inextensible cloth model We use four different meshings of the cloth: two with triangles and two with quadrilaterals (see Figure 2.4), but we keep the physical parameters fixed. The visual results of the experiment can be seen in Figure 2.5. The four textiles fold diagonally without any locking artifact. Figure 2.5: Four different meshes used for the experiment: the red ones are made of triangles and the blue ones of quadrilaterals. From left to right the number of nodes is: 472,941,529,900. On the bottom of each cloth, we display the euclidean position of the free corner. nodes /element 472/△941/△529/□900/□ 472/△0 2.45 cm 3.16 cm 1.41 cm 941/△-0 1.41 cm 2.83 cm 529/□- - 0 3.74 cm 900/□- - - 0 Table 2.2: Distances (cm) between the bottom corner of the cloths (of dimensions 1m×1m). In Table 2.2we compute the euclidean distance of the bottom corner of every textile with respect to each other. Note how, doubling the number of nodes, or changing from irregular triangles to regular quadrilaterals, does not alter the result within a margin of error of a few centimeters (recall that the sheets are 1m×1m). 2.7.2Cusick’s test This subsection will demonstrate the stability of our model under refinements of the mesh. We will explain how to simulate Cusick’s test [56] using our simulator. Cusick’s test consists in letting a circular cloth of radius 15cm drape on top of an also circular table of radius 9cm (see Figure 2.6) with their centers aligned. 2.7 evaluation and results 31 Figure 2.6: Photo of a circular cloth draping on top of a circular table. Image taken from [42]. Then, one computes the area of the plane projection of the draped cloth C and divides it by the area of the flat ring A of width 6cm (see Figure 2.7). Figure 2.7: After the cloth has draped on top of the table, the drape coefficient is calculated as the ratio between the area of the plane projection of the cloth C , divided by the area of the flat ring A of width 6cm. Image adapted from [42]. This is called the drape coefficient of the cloth: DC % =100 ×C A. (2.22) It is known that this coefficient depends heavily on the stiffness of the cloth [56]. This can be understood intuitively: a completely rigid sheet of metal would not bend and would have DC equal to 100% , whereas a really light material (e.g. silk) has a DC pretty close to 0 . The measurement of this coefficient is not trivial: it usually requires a dedicated machine (Cusick’s Drape Tester [65]); and variation of the measured coefficient for the same cloth is common, so several specimens of the same textile have to be employed and then their DC’s averaged [42]. Naturally, these complications have caused the creation of alternative measuring methods using computer simulations (e.g. [38]). In order to simulate Cusick’s test with our model, we quadrangulate the ring A of Figure 2.7, we fix in space its inner boundary nodes and let gravity act (see Figure 2.8). In order to study the dependency of the DC with respect to mesh resolution, in Figure 2.9we plot the DC for three stiffness values κ (low: κ0 , medium: 5κ0 and high: 10κ0 ; where κ0=0.005 ) and 24 different meshes of the ring A ranging from 250 nodes to 1300 . The three mean values for the computed DCs are 32 inextensible cloth model Figure 2.8: Simulation of Cusick’s test with the inextensible model and a mesh of 768 nodes. From left to right we vary the stiffness of the cloth and get respectively DCs of 26.1%, 54.9%, 77.2% . They correspond to peach-skin polyester (low stiffness), imitation wool (medium) and tussore cotton (high). 23.2%, 52.4%, 77.3% and they approximately correspond to peach-skin polyester (low stiffness), imitation wool (medium stiffness) and tussore cotton (high stiffness) (see the Appendix of [42]: IDs: 38, 37, 22 ). Figure 2.9shows that the computation of the DC with our model is very stable (especially from 700 nodes on), and hence the model has a robust behavior with respect to mesh resolution. As we already said, this is very relevant for robotic applications, where the use of coarse meshes becomes a necessity because of performance constraints. Figure 2.9: Computation of the drape coefficient (DC) for 3stiffness values (low, medium and high) and 24 different meshes of cloth from 250 to 1300 nodes. The three mean values for the computed DCs are 23.2%, 52.4%, 77.3%. 2.7.3Tank-top simulation With this experiment, we aim to show the ability of our model to simulate complex topologies (e.g. a tank-top shirt, see Figure 2.10 and 2.7 evaluation and results 33 https://youtu.be/DvSWEXdw6Bo ) and study different area errors. We track the variation of area of the garments through time. We think area is a good measure since it is a physical property of cloth and not a mesh-dependent metric (e.g. stretching of the edges of the elements). We compute the area of each element of the discretized cloth using its local parametrization (2.4). Figure 2.10: Simulation with the inextensible model of the shaking of a triangulated tank-top. At each time instant ( t=1.1s, 1.75s, 2.5s, 3.5s ) we plot the area error with sign (2.26) of each individual triangle. We use three area errors: 1.Total area error: et(t) = |A(t)−A0| A0 , (2.23) where A0 is the total area of the garment at time t=0 and A(t) is its area at time t>0. 2.Mean element error: em(t) = 1 ne∑ e |ae(t)−ae(0)| ae(0), (2.24) where ae(t)is the area of element Ωeat time t≥0. 3.Dispersion: da(t) = em(t) + 2sVare|ae(t)−ae(0)| ae(0). (2.25) We simulate the shaking with a hanger of a meshed tank-top with 1676 triangles (see Figure 2.10) during 4.5s. There, we also plot the area error with sign of each individual triangle for each time instant: e(t) = ae(t)−ae(0) ae(0). (2.26) It is interesting to notice the concentration of error near the boundaries of the cloth and at the points where the tank-top makes contact 40 cloth collisions takes advantage of the way we integrate numerically the equations of motion. For the time being, assume we have both φn and φn+1 (and their velocities) available to make computations. 3.3.1Detection of self-collisions In general, we assume that the cloth is triangulated (in case of a quadrangulation we can always divide the quads in two); then in case of collision, there are only two stable (i.e. detectable) possibilities: an edge-edge collision and a node-face collision. In these two cases, we have four nodes involved which at some instant of time belong to the same plane (see Figure 3.1). We must then only check if two co-planar edges cross or if a point belongs to a triangle. These two problems are readily solved using barycentric coordinates. Figure 3.1: If the moving triangles were not intersecting before, then there exists a time in which the edges ab and a′b′were co-planar. Now we describe in more detail the process: in order to save computational time, we only check if a collision has happened for pairs of edges (or nodes and faces) that at time tn are sufficiently close (and not for every pair, which would result in combinatorial explosion). To obtain this list of sufficiently close (up to some tolerance) pairs, there are several methods: one of the most widely used is the hierarchical method [92], where the mesh is divided into large regions, which are in turn also divided in smaller regions, up to a certain number of times (the number of hierarchies). Then one only checks if two pairs are sufficiently close when all of the larger regions containing these are close enough (e.g. by computing the distance between their center of masses), otherwise, they are discarded. Since our meshes are fairly coarse we avoid hierarchies (which can be cumbersome to implement) and compare the center of masses of our pairs in order to detect those that are candidates for collision. Next, denoting by x1,x2,x3,x4 the position of the four candidate nodes at time tn and by v1,v2,v3,v4 their velocities, we must only check if 3.3 self-collisions 41 det(˜ x1+t·˜ v1,˜ x2+t·˜ v2,˜ x3+t·˜ v3) = 0 where ˜ xi=xi−x4 and ˜ vi=vi−v4 , since φn+1=φn+dt ·˙φn+1 . This is a cubic equation a3t3+a2t2+a1t+a0=0 in t, with coefficients:                a3=det(˜ v1,˜ v2,˜ v3); a2=det(˜ x1,˜ v2,˜ v3) + det(˜ v1,˜ x2,˜ v3) + det(˜ v1,˜ v2,˜ x3); a1=det(˜ x1,˜ x2,˜ v3) + det(˜ x1,˜ v2,˜ x3) + det(˜ v1,˜ x2,˜ x3); a0=det(˜ x1,˜ x2,˜ x3); (3.3) When t≪dt is small, the solution of the previous equation can be approximated linearly by −a0 a1 . In any case, if there is a root for some tc∈[0, dt] , we must then do two different calculations with the four co-planar points yi=xi+tc·vi, namely: 1. Edge-edge case: say the endpoints of the first edge are y1,y2 and of the second y3,y4 , then we compute the real numbers α,β where y1+α(y2−y1) = y3+β(y4−y3) (the intersection of the two lines defined by the segments) and if they belong to the interval [0, 1]a collision has occurred. 2. Node-face case: say the node is y4 and the other 3points are the corners of the triangle, then we compute the real numbers u,v,w such that y4=uy1+vy2+wy3 and if they are in the interval [0, 1]a collision has occurred. 3.3.2Constraint definition for self-collisions We now describe the computation of the self-collision constraint Hk . It will be linear in φ and naturally have slightly different forms depending on our two cases: 1. Edge-edge case: Hk(φ):=⟨πα(x1,x2)−πβ(x3,x4),ν⟩ ≥ 0, where xi are the four endpoints of the two edges, πα(x1,x2) = (1−α)x1+αx2 and πβ(x3,x4) = (1−β)x3+βx4 are the closest points between the two segments and ν is the normal vector to both edges. In general, the values ν,α,β vary with time. We will nevertheless assume that they are constant during the time-step, and compute them with the positions of the segments given by φn+1. The normal vector νis oriented such that Hk(φn)≥0. 2. Node-face case: Hk(φ):=⟨x4−π(x1,x2,x3),ν⟩ ≥ 0, 42 cloth collisions where x4 is the node, xi are the 3corners of the triangle, again π(x1,x2,x3) = ux1+vx2+wx3 is the closest point inside the face to the node and ν is the normal vector to the triangle. In general, the values ν,u,v,w vary with time. We will again assume that they are constant in time, and compute them with the positions given by φn+1 . The normal vector ν is oriented such that Hk(φn)≥0. Remark 3.3.1.Notice that: 1. By construction Hk(φn+1)<0. 2. The constraint Hk is an approximation of the signed distance between the pairs edge-edge and node-face (only an approximation since νand the barycentric coefficients are fixed in time). 3. Since in practice cloth has thickness, say τ0 , the constraint we actually must impose is Hk(φ)≥τ0. 3.3.3Proximity constraints and cloth thickness Adding the constraints we have just defined is enough to correct all present self-collisions. Nevertheless, there are two main drawbacks: 1. Efficiency: most cloth collisions can be avoided before they happen by adding preventive constraints. 2. Vibrations: since we are assuming that the cloth has a thickness τ0>0 , when we integrate the system again and go from Hk(φn)<0 to Hk(φn+1)≥τ0 , the change between the position of the nodes can be too large, and since our cloth is inextensible, this could create unwanted oscillations. In order to avoid these two problems, we apply the detection procedure previously explained in 3.3.1with one small difference: during the detection phase we move the pairs (edge-edge or face-node) closer, using their normal vectors and taking into account the thickness of the cloth, so that pairs that are too close and/or are approaching each other, are kept at a minimum distance of τ0 before they actually cross. Since the restrictions we are considering are inequalities, we can add these to the system because they only affect the dynamics of cloth in case the constraint will actually get violated. In symbols, this means that we compute the coefficients (3.3) of the third degree polynomial using the altered positions given by ˆ xi=xi±ωτ0ν , where ν is the unit normal vector (the cross product for the edge-edge case and the normal to the triangle for node-face case), ω≈0.5 is what we will call a proximity parameter and the sign ±is chosen so that the pairs approach each other. Afterwards the response constraint Hk is calculated as usual (i.e. the normals and the barycentric coordinates) with the unaltered positions xigiven by φn+1. 3.4 numerical integration of the system 43 Remark 3.3.2.It is usually enough to use the positions ˆ xi=xi±ωτ0ν , where ω≈0.5 to detect all self-collisions, nevertheless some can sometimes be missed because the nodes have moved too much. In that case, we enter an iterative process reducing gradually the value of ω until all are resolved. We will explain this in more detail in Section 3.4.1. Definition 5.(Self-collision constraints). We will denote by C=Collisionsωφn→dt φn+1 the set of self-collisions constraints that must be imposed from the state φnto the state φn+1with proximity parameter 1 2>ω≥0. 3.4 numerical integration of the system The friction force and the contact constraints introduced in the previous section are in general highly non-linear and stiff and thus must be integrated implicitly (like the inextensibility constraints, see Section 2.6.1). To integrate the system numerically from time tn to tn+1 (i.e. to advance the simulation from φn to φn+1 ), we perform as before an iterative process φj+1=φj+∆φj+1 where the initial point is the unconstrained step φ0=φn+1 0(φn, ˙φn) given by an implicit Euler scheme. Also, we write: H(φj+1) = H(φj+∆φj+1)≃H(φj) + ∇H(φj)∆φj+1, and similarly C(φj+1) = C(φj+∆φj+1)≃C(φj) + ∇C(φj)∆φj+1, and then solve iteratively the following sequence of quadratic programs with linear equality and inequality constraints:          min∆φj+1 1 2∆φ⊺ j+1·M·∆φj+1−∆φ⊺ j+1·fµ(˙φj) C(φj) + ∇C(φj)∆φj+1=0, H(φj) + ∇H(φj)∆φj+1≥0, (3.4) where 1. ˙φj+1=φj+1−φn dt is an approximation of ˙φn+1, 2.fµ(˙φj) = −µV(˙φj)⊺∆βjis the friction force at iteration j, 3.V(˙φj)are the relative unit tangent velocities, 4.(∆βj)i=||∇Hi(φj)⊺(∆γj)i|| is the magnitude of the contact forces at iteration j, 44 cloth collisions 5. and ∆γj≥0 are the multipliers associated to the contact constraints. We iterate until max |C(φj)|<ϵ0, min H(φj)≥ −ϵ1, max |∆φj|<ϵ2(3.5) for some tolerances ϵ0,ϵ1,ϵ2>0 . This third condition ensures that the friction force has stabilized. Note that the critical points of the previous quadratic problems (3.4) are:                        M·∆φj+1=−∇C(φj)⊺∆λj+1+∇H(φj)⊺∆γj+1−µV(˙φj)⊺∆βj, C(φj) + ∇C(φj)∆φj+1=0, H(φj) + ∇H(φj)∆φj+1≥0, ∆γj+1≥0, ∆γ⊺ j+1·hH(φj) + ∇H(φj)∆φj+1i=0, (∆βj)i=||∇Hi(φj)⊺(∆γj)i||. (3.6) Remark 3.4.1.In order to integrate friction force we have made the approximation fµ(˙φj+1)≃fµ(˙φj) . That is, we have dropped the gradient we would normally have with a first-order approximation (this is what we also do with the gradient of the constraint forces in the first equation of (3.6)). 3.4.1Addition of self-collision constraints Instead of checking and generating all self-collision constraints only with the states C=Collisionsωφn→dt φn+1 , we take advantage of the fact that we perform an iteration process. We now explain how we introduce self-collisions into the sequence of problems (3.4) for every step. For every iteration j we check for self-collisions (see Section 3.3.1) taking into account the thickness of the cloth (Section 3.3.3) between the states φn and φj and generate the corresponding constraints (Section 3.3.2). In symbols this means that all the constraints Cj=Collisionsωφn→dt φj for j≥0 and ω≈0.5 are added to the system. In the rare case that the same collision is found in two different iterations we only keep the constraint defined by the later iteration. Then, when we find a state φj+1 that satisfies the stopping criteria (3.5), we check for self-collisions with ω=0 , and in case no self-collision is detected, we put φn+1=φj+1 . Otherwise, we repeat the whole iteration process with a smaller value of ω (see Remark 3.3.2). 3.5 efficient solution of the quadratic problems 45 3.5 efficient solution of the quadratic problems Definition 6(Active constraint).In a constrained optimization problem (such as (3.4)), we say that an inequality constraint g(x)≥0 is active at a feasible point y if g(y) = 0 . Note that all equality (in our case inextensibility) constraints are always active. In order to solve the sequence of problems (3.4) we could employ any quadratic problem solver, but we would not be taking advantage of the structure of our problem. That is, if in one of the iterations j one of the contact constraints Hi is active (see the previous definition), then it is likely that it will be active again at the next iteration. Physically, this means that nodes of the cloth that are in contact with an obstacle (or among themselves) at some iteration, are likely to remain in contact. This suggests the use of active-set-methods [86] to solve the quadratic problems. We will develop a novel active-set algorithm in the following pages. Although we could use one of the many existing ones, they always require that one begins with a feasible (albeit not optimal) solution to the problem. Our method will not have this requirement. The main idea of active set methods is to find the active set of constraints at the solution, because then, once known, the program can be solved by ignoring inactive constraints, and assuming that all active inequality constraints are equality constraints. Recall that solving quadratic problems with equality constraints is very cheap and can be done by solving a linear system (see [44]). This will be precisely what we will do for every iteration of the sequence (3.4). In order to find the active set, one splits the constraints in two sets: The working set, W : these are the constraints believed to be active ( g=0 ) and therefore are imposed as equality constraints when one solves the optimization problem. This can be initialized as the set consisting only of equality constraints. The observation set, O : these are the constraints believed to be inactive ( g>0 ) and therefore are not imposed as equality constraints. Since they are not included in the problem one must be careful that they do not become violated. Then one proceeds as follows: 1) solve the equality problem defined by the working set; 2) compute the Lagrange multipliers of the working set for the inequality constraints; 3) send a subset of the constraints with negative Lagrange multipliers to the observation set; 4) if all multipliers are positive, check if all constraints in the observation set remain feasible; 46 cloth collisions 5) send a subset of the infeasible constraints to the working set; 6) repeat. Then, if at some iteration we have found an increment ∆φj+1 such that all contact constraints in the working set have positive Lagrange multipliers ∆γj+1≥0 (see Equation (3.6)) and all constraints in the observation set are not violated, we have found the active set (see [86]) and we can make the update φj+1=φj+∆φj+1 . The following proposition ensures that we do not enter in a never-ending cycle: Proposition 1(Entry and exit of constraints).Given the system (with unknowns ∆φ)    M∆φ=∇H(φ)⊺∆γ, H(φ) + ∇H(φ)∆φ=0, (3.7) and the system (with unknowns ∆˜φ)    M∆˜φ=∇H−k(φ)⊺∆˜γ, H−k(φ) + ∇H−k(φ)∆˜φ=0, (3.8) where we have removed the constraint Hk(φ) + ∇Hk(φ)∆φ=0 from the first system; then it holds that ∆γk·(Hk(φ) + ∇Hk(φ)∆˜φ)≤0, (3.9) where ∆γkare the Lagrange multipliers of the removed constraint k. Proof. Subtracting the first two equations of the systems, we get: M(∆˜φ−∆φ) = ∇H−k(φ)⊺(∆˜γ−∆γ−k)− ∇Hk(φ)⊺∆γk. Then, multiplying both sides by (∆˜φ−∆φ)⊺, we deduce that 0≤(∆˜φ−∆φ)⊺·M·(∆˜φ−∆φ) = 0−(∆˜φ−∆φ)⊺· ∇Hk(φ)⊺∆γk, since (∆˜φ−∆φ)⊺· ∇H−k(φ)⊺=H−k(φ)⊺−H−k(φ)⊺=0. Finally, using that ∇Hk(φ)∆φ=−Hk(φ), and rearranging terms we get 0≤ −∆˜φ⊺∇Hk(φ)⊺∆γk−Hk(φ)∆γk. From here (3.9) follows easily. Corollary 1.If a constraint in the working set has a negative Lagrange multiplier, when it is taken out of the system and put in the observation set, it becomes feasible. Conversely, when a constraint in the observation set is infeasible and we sent it to the active set, its associated Lagrange multiplier is positive. 3.5 efficient solution of the quadratic problems 47 Remark 3.5.1.The heuristic that is usually followed to decide which constraint to remove or to add is: delete from the working set the constraint with the most negative Lagrange multiplier and add to the working set the constraint from the observational set that is being most violated (the most negative one). To finish this section we study the case of linearly dependent constraints. This is relevant since in general, we do not want to introduce linearly dependent constraints into the system because they give raise to (near) singular matrices. Lemma 1.If a constraint G in the observation set can be written as a linear combination of constraints of the working set, i.e. G(φ) = ∑αiHi(φ) , then the linearized constraint is feasible G(φ) + ∇G(φ)∆φ≥ 0. Proof. G(φ) + ∇G(φ)∆φ=G(φ) + ∑αi∇Hi(φ)∆φ=G(φ)−∑αiHi(φ) = 0. Remark 3.5.2.The previous lemma ensures that in general, we do not send linearly dependent constraints from the observation set to the working set. Nevertheless, it is possible to have a degenerate case where the constraints are not linearly dependent but their gradients are. In symbols, this would mean that a constraint in the observation set satisfies G(φ) + ∇G(φ)∆φ≤0 and moreover ∇G(φ) = ∑αi∇Hi(φ) . What we do then is to introduce G in the working set while removing the Hi with the greatest αi=0 in absolute value. This new working set is linearly independent (otherwise it would contradict the assumption that the original working set without Gwas linearly independent) and the process can continue. 3.5.1Factorization of the matrix system Every time that a constraint goes from the working set to the observation set (or viceversa), i.e. when the index sets W and O are updated, the Lagrange multipliers must be recomputed, i.e. a linear system must be solved in order to find the solution of (3.10).          M·∆φj+1=−∇C(φj)⊺∆λj+1+∇H(φj)⊺∆γj+1−µV(˙φj)⊺∆βj, C(φj) + ∇C(φj)∆φj+1=0, Hi(φj) + ∇Hi(φj)∆φj+1=0 for i∈ W. (3.10) In order to ease readability we will include inextensibility constraints and the contact constraints of the working set in only one 48 cloth collisions function denoted by G⊺= [C⊺,H⊺] . Now, since in general only one constraint will be entering or exiting at the time, the linear systems that we have to solve are almost identical with the exception of a few rows and columns. That is why, the use of factorizations becomes an important tool to achieve efficiency. The linear system we need to solve to find the multipliers is: ∇G(φj)M−1∇G(φj)⊺∆ζj+1=−G(φj)− ∇G(φj)M−1fµ(φj), where ∆ζ⊺ j+1= [∆λ⊺ j+1,∆γ⊺ j+1] . Hence the system matrix (let us call it A ) is positive definite (since M is positive definite because it is the mass matrix); therefore we can use Cholesky decomposition [45], provided our constraints are linearly independent (see again Lemma 1 and Remark 3.5.2). This is a factorization of the form QAQ⊺=LL⊺, (3.11) where L is an invertible lower triangular matrix and Q is an (orthonormal) permutation matrix (this is done in order to take advantage of sparsity patterns). Solving the linear system in this way reduces to a couple of triangular substitutions and two matrix multiplications: A−1= (Q⊺LL⊺Q)−1=Q⊺L−1(L⊺)−1Q, since L and L⊺ are lower and upper triangular, and we can use forward and backward substitutions (no matrix is actually inverted). 3.5.2Updates of the Cholesky decomposition Every time a constraint enters or exits the working set, L (and Q ) can be efficiently updated without recomputing the factorization from scratch (see, e.g. [28,102]). In the following, we describe two simple methods (not necessarily the most efficient) to perform such a task. When introducing constraints, we will do so, only one at a time, whereas to remove them we will usually delete all constraints with negative Lagrange multipliers. Assume first that we want to add a single contact constraint Hk to the working set. Then the updated matrices are: L+="L 0  l⊺λ#,Q+="Q 0  0 1#, (3.12) where lis found by solving L· l=Q· ∇H·M−1·(∇Hk)⊺and λ=q∇Hk·M−1·(∇Hk)⊺− l⊺· l This can be easily seen by writing the expanded equations for: 3.5 efficient solution of the quadratic problems 49 Q+A+Q⊺ +=L+L⊺ +, (3.13) where A+="∇G ∇Hk#·M−1·"∇G ∇Hk#⊺ . (3.14) Remark 3.5.3.The case when λ≤0 occurs when we try to introduce a linearly dependent constraint into the system. We use this as a test to check if we must interchange constraints between the working and observation set as explained in Remark 3.5.2. On the other hand, updating Lwhen taking out an equation is more complicated since when we remove a row (or a column) from L it loses its triangular form. There are several methods to accomplish this update without recomputing everything from scratch, among them: using low-rank downdates [45], using artificial multipliers to set to zero the solution’s coordinates corresponding to exiting constraints or employing Givens’ rotations (also known as Householder transforms, see [102]). Among these three methods, we have found that the most practical and fastest for vectorized languages such as Python or MATLAB is the artificial multipliers method, which we now describe. Assume we have the system: Ax =b, where every coordinate k of x is the Lagrange multiplier corresponding to a constraint Gk . Then, if we want to remove the constraints i1, . . . , in from the system, we set their values to zero and hence solve the system    Ax +S⊺y=b, Sx = 0, (3.15) where Sk·x=xik . This is the case since equations (3.15) are the critical points of:    minz1 2z⊺Az −z⊺b Sz = 0. Remark 3.5.4.As already mentioned, recomputing from scratch the full Cholesky factorization is usually not desirable, especially in cases with very fine meshes. Nevertheless, when the number of variables is not too high, it can be competitive in the case of the removal of constraints, since we can delete more than one condition a time without compromising linear independence (as opposed to the case of adding constraints, where adding more than one at the time can be problematic). With this method, the heuristic used is that all constraints with negative Lagrange multipliers are removed from the working set. 56 cloth collisions The values ν,u,v,w are (like in the case for self-collisions) constant in time, and are computed with the positions given by φj . The normal vector νis oriented such that Hi(φn)≥0. Remark 3.6.1.As with self-collisions we consider cloth’s thickness in practice by imposing Hk(φ)≥τ0>0 . Moreover, as before this thickness is taken into account in the detection process (see Section 3.3.3). In Figure 3.4we can observe the result of the simulation from three different viewpoints. The cloth lies stably on top of the cusps without any noticeable artifact. In Chapter 4, we will encounter again sharp objects when we simulate the hitting of a piece of cloth with a thin stick. 3.6.4Folding sequence of short pants: https://youtu.be/2gdnjUICb0g In this final experiment, we simulate the dynamical folding of a pair of shorts (the same ones of Figure 2.2) on top of a table. In order to do so, we control two nodes at the top of the shorts. In this experiment, we also allow the garment to shear by introducing the force described in Section 2.4.6. Figure 3.5: Trajectory of the controlled nodes for the dynamic folding of shorts. The first part of the motion is performed fast enough so that the shorts have sufficient momentum to lay partially flat on top of the table after lowering them. Finally, the fold is completed by dropping the top two corners on top of the leg loops. The followed trajectory of the controlled nodes can be seen in Figure 3.5. 3.6 evaluation and results 57 Figure 3.6: Simulated sequence of the dynamical folding of a pair of shorts. The first part of the motion (frames two and three) is performed fast enough so that the shorts lay partially flat on top of the table after a lowering phase (frame four). The final fold is completed by dropping the top two corners on top of the leg loops (frames five and six). In Figure 3.6we depict six frames of the simulation. Notice how crucial is the well-functioning of the self-collisions algorithm for a realist outlook of the whole folding sequence. 4 EXPERIMENTAL VALIDATION The definitive test for a cloth’s model is its comparison to reality. Textile engineering models are focused on such comparison, to the point of developing specialized testing equipment. But the object of their study are local properties of cloth, such as elasticity parameters, which are tested in static scenarios (e.g. [24,56,79,117]). Other, more recent lines of research [96,97] focus on estimating friction coefficients using non-intrusive video images. To the knowledge of the authors, none of the models previously mentioned has been able to compare its results with the motion of cloth dynamically. In this chapter, we seek to study how faithfully our inextensible model can reproduce recordings of real textiles under different circumstances. This chapter is divided in 3parts: WAM-shaking experiments (Section 4.1): in this first round of experiments we record the motion of four (size A3) textiles with a depth camera. The fabrics are shaken by a robotic WAM arm using two sets of different amplitudes and frequencies. The goal is to assess with how much accuracy our model is capable of reproducing the dynamics of the textiles. Aerodynamics study (Section 4.2): in this second round of experiments we record the motion of eight (four size A3and four size A2) textiles with a Motion Capture System. The fabrics are shaken and twisted by a human at two different speeds. We carry out two repetitions of each motion and therefore have 64 different recordings of about 15 seconds. The goal is to study how the speed and size of the textiles affect their motion and afterwards develop a predictive and a priori formula for the value of the cloth’s physical parameters. Collision’s validation (Section 4.3): in this final set of experiments we use again motion capture, and record the collision of four (size A2) textiles. In one of the scenarios, the fabrics are laid dynamically on top of a table in a putting-a-tablecloth fashion. In the other, they are hit by a long stick four times at various places and with different strengths. The goal is to assess the accuracy of the collision and friction model previously developed and put to use the predictive formulas found before. 59 60 experimental validation Cloth’s materials and sizes For the experiments in this chapter, we employ seven cloth materials (see Figure 4.1) and two different sizes: A3(0.297 x0.420 m with area 0.1247 m2 ) and A2(0.42 x0.594 m with area 0.2495 m2 ). Before performing the experiments they were ironed to remove all considerations of plasticity from the validation process. Not all the textiles are used in all experiments, in every scenario, we will specify what materials and sizes are used. In Table 4.1we can see the density of all the fabrics and some typical examples of garments made from them. Figure 4.1: In the picture we can see all the fabrics (size A3) used in the experiments. From left to right we have: paper, polyester, light cotton, felt, wool, denim and stiff cotton. fabric density (kg ·m−2)sizes examples Paper 0.0802 A3Polyester 0.1042 A3, A2Silk-like. Light-cotton 0.1604 A3Dressing shirt. Felt 0.1764 A3Wool 0.1804 A3, A2Formal suit. Denim 0.3046 A3, A2Jeans. Stiff-cotton 0.3046 A3, A2Sack. Table 4.1: Density, sizes and examples of all the materials used in all the experiments. 4.1 wam-arm experiments 61 4.1 wam-arm experiments A Barret robotic arm together with a depth camera are employed to record the real motion of garments being shaken at different velocities (see Figure 4.2). We implement oscillatory motions: a forward and backward shaking (see Figure 4.5) at two different speeds. Afterwards, the resulting point-cloud recordings are de-noised, interpolated and meshed, so that we end up with the spatial trajectory of the vertices of a polyhedron. In order to keep this first validation manageable, we focus only on four different size A3textiles: light cotton, wool, felt and paper. We use a trivial rectangular topology in order to avoid occlusions in the point-cloud when using the depth camera to record the motions. Finally, we find the physical parameters of the model that best approximate the motion of the polyhedron using a least squares approximation and a minimization algorithm. We study the evolution of absolute errors and dispersion measures (i.e. standard deviations) of the simulated garments versus the real ones for all textiles and motions. In the following sections, we give more specific details of the whole recording and comparison process. Figure 4.2: Experimental setup for the recording of the motion of the real textiles. On the left, the depth camera. On the right, the robotic arm with an affixed hanger. 4.1.1Camera data As mentioned before, comparison with reality is performed by recording in the laboratory the motion of a piece of cloth when subjected to exactly the same movements of the robotic arm controlling it as specified in the simulations with our model. The real motion is captured 62 experimental validation by a depth camera Kinect XB360-607 (see Figure 4.3). The garment is a rectangular piece of cloth, with its two upper corners fixed to a rigid hanger, which is affixed itself to a robotic arm Barrett WAM. To keep the garment completely in view of the camera, the robot moves it along the depth axis of the camera (axis x in our reference) following a curve of equation: (x(t),y(t),z(t)) = (Acos(2πft) + c,y0,z0). (4.1) The two corners of the garment fixed to the robot arm follow the same oscillation, with different but constant values y1,y2 instead of y0 . Figure 4.3: Point-cloud obtained using a depth camera (bottom) and RGB image (top). Note that the camera only gets the depth of the objects visible to it, all the rest (e.g. the robot behind the cloth) is occluded. Two motions are recorded: the slow one has A=0.15 m and f= 0.3 Hz, while the fast one has A=0.075 m and f=0.6 Hz. Surrounding objects are filtered using planes for each recorded frame, and outlying points are detected on the basis of distance and removed. Only garments with a homogeneous color are tested, so we can remove further noise through the use of a color filter. After the point-cloud has been thus filtered, it is meshed using the algorithm described in [55] (see Figure 4.4). The model that we are validating is continuous, producing 4.1 wam-arm experiments 63 simulations that are very stable under remeshing. Because of this, it suffices to select as a baseline for all the experiments a coarse 9×9 mesh. Figure 4.4: Quadrilateral meshing (left) of one of the frames of the filtered and de-noised point-cloud (right). Each quadrilateral is divided into triangles for plotting purposes. The robotic arm follows the oscillation of Equation (4.1) for 10 seconds. The camera starts recording 1 second before the start, and finishes 3 seconds after the end of the robotic motion, for a total of 14 seconds of recorded garment motion. The original recording has its frames at irregular time steps. We replace them through linear interpolation to obtain the meshed cloth’s position each dt =0.01 seconds. 4.1.2Parameter fitting We will be adjusting only 2parameters: the damping parameter α and the (virtual) gravitational mass δ . The other parameters in Table 2.1 are set to 0 except for ρ=1 , since α and δ are the most relevant for the motions performed in this experimental validation ( β and κ affect mostly the local deformation of the cloth, see [7]). In Section 4.1.1we explained how to obtain a sequence of positions of the nodes of the real cloth {ϕ0,ϕ1, . . . , ϕm} . If we integrate numerically Equation ( 2.19) using the same trajectories (4.1) for the two upper corners, we get a sequence {φ0(δ,α),φ1(δ,α), . . . ,φm(δ,α)} of positions of the nodes of the simulated cloth (where φ0=ϕ0 ) for each value of the two parameters. Hence, a natural error metric to minimize is: L(δ,α) = [m 2] ∑ i=1 ||φi(δ,α)−ϕi||2 M, (4.2) where || · ||M is the L2 norm with respect to the matrix M (i.e. ||x||2 M= x⊺·M·x ), and we only use the first half (7seconds) of the recorded frames to perform the fitting. We call this the training window. Finally, we minimize L using a derivative-free algorithm (the Nelder-Mead Simplex Method) and find the optimal values of δ,α . The resulting parameters are shown in Table 4.2. 64 experimental validation Remark 4.1.1.Notice that for these experiments we have taken ρ=1 and hence the optimal values of δ,α will not be comparable to those found in the other next experiments when we use the real values of ρ of Table 4.1. We will analyze the evolution of the time-dependent absolute error: ei(δ,α) = q||φi(δ,α)−ϕi||2 M. (4.3) material δslow δf ast αslow αf ast ¯ eslow ¯ ef ast Paper 0.32 0.37 1.40 2.49 0.45 cm 0.41 cm Light-cotton 0.52 0.52 1.29 2.69 0.41 cm 0.31 cm Felt 0.46 0.61 1.05 2.65 0.39 cm 0.30 cm Wool 0.47 0.50 1.19 2.52 0.37 cm 0.34 cm Table 4.2: Estimated parameters and mean absolute errors (cm) for the slow and fast oscillatory motions performed by the WAM robot. In the last two columns of Table 4.2we show the average values of ei over the testing window (i.e. the last 7seconds of the motion). Moreover, for plotting purposes we use two times the standard deviation of the error at each node jas a dispersion measure: di(δ,α) = ei(δ,α) + 2rVarj∈Nodes(S)||φi j(δ,α)−ϕi j||R3. (4.4) Remark 4.1.2.Notice that the variance is taken along the nodes of the surface and not on time, hence we are measuring a spatial standard deviation. Figure 4.5: Comparison at four time instants of the recorded fast motion of paper (right) versus its simulation with the inextensible model (left) with δ=0.37 and α=2.49 . The mean absolute error is 0.41 cm. 4.1 wam-arm experiments 65 4.1.3Sensitivity analysis In order to show how robust our model is when away from the optimal values found in Table 4.2we perform a sensitivity analysis. We carry out 100 different simulations with different values of the parameters α≥0 and δ>0 (with an upper limit of two times their optimal value), and compute the mean of the absolute error (4.3) on the testing window formed by the last 7seconds of the movement. The results for the fast motion of paper (see Figure 4.5) are shown in Figure 4.6. Figure 4.6: Mean of the absolute error (4.3) for the fast motion of paper using different values of the parameters of the model (damping α and virtual mass δ ). In red we highlight the error with the optimal value of the parameters. First of all, note that for the discrete set of values considered, the optimal values found are in fact a global minimum. It is interesting to note that the model is very robust with respect to the virtual mass δ , and in fact even when δ=1 (no aerodynamics), there is a value of α=2.99 that gives a very low overall error (0.61 cm). Other interesting limit cases are α=0 (no damping) and α=5 (over-damped), but they give higher errors. Nevertheless, we do not consider δ=0 (no gravity), since it is very unrealistic, and that is confirmed in the picture since the largest errors occur around δ=0.05. 72 experimental validation 4.2.1Movements and textiles For this second set of experiments, we employ the following materials: polyester, wool, denim and stiff-cotton. Both A2and A3sizes are used. For the A2textiles we use 20 reflective markers, whereas for the A3 ones 12 are used. In both cases, the makers are placed equidistantly in order to obtain a faithful representation of the dynamics of the fabrics. In contrast to the first experiment with the WAM robot, this time the motions are performed by a human. This introduces way more uncertainty than before, since every movement has its own unique variabilities. Therefore every motion was recorded twice on different days: repetition I with a special hanger (see Figure 4.11) and repetition II with bare hands (see Figure 4.12). The motions are: Figure 4.13: Shaking motion sequence (left to right): the cloth is shaken back and forwards. 1. Shaking: this is similar to the motion performed with WAM robot, the cloth is shaken back and forwards (see Figure 4.13). Figure 4.14: Twisting motion sequence (left to right): the cloth is rotated with respect to the z-axis back and forth several times. 2. Twisting: this is a new motion where the cloth is rotated multiple times (approximately 30 degrees) with respect to the z -axis (see Figure 4.14). Each motion lasts approximately 15 seconds (with a frame every dt =0.01 seconds) and is performed at two different speeds: slow and fast. In Table 4.4we can see some values of average speeds ( m·s−1 ) comparing fast and slow motions for the twisting movement of the A2textiles. 4.2 aerodynamics study 73 material slow i fast i slow ii fast ii Polyester 0.081 0.250 0.090 0.247 Wool 0.093 0.228 0.085 0.256 Denim 0.092 0.308 0.090 0.292 Stiff-cotton 0.091 0.218 0.107 0.246 Table 4.4: Average velocities ( m·s−1 ) for the twisting motion of the A2textiles. We display the speeds for the first repetition (with a hanger) and for the second (with bare hands). We can of course observe some variability but overall the speeds are maintained pretty consistently. Remark 4.2.2.Notice that we have two motions at two speeds for four textiles with two different sizes repeated two times, which makes a total of 64 recordings. 4.2.2Parameter fitting As before we will be adjusting only 2parameters: the damping parameter α and the (virtual) gravitational mass δ . The other parameters in Table 2.1are set to 0 except for the density ρ , which this time is set to its corresponding value in Table 4.1. We denote the sequence of positions (this time given by following the reflective markers) of the recorded fabric’s nodes by {ϕ0,ϕ1, . . . , ϕm} . Again, we integrate numerically Equation ( 2.19) using the recorded trajectories of the two upper corners (i.e. the two markers placed at the top corners), and get a simulated sequence {φ0(δ,α),φ1(δ,α), . . . ,φm(δ,α)} (where φ0=ϕ0 ) for each value of the two parameters. For the simulations we consider a refinement of the initial meshes given by the markers: for the A3case we employ a 5×7 meshing and for the A2case a 7×9 one. As metrics we use again the absolute error (4.3) (only at the recorded nodes) and the following time-dependent (on i ) spatial (on j ) standard deviation: si(δ,α) = rVarj∈Nodes(S)||φi j(δ,α)−ϕi j||R3. (4.5) Remark 4.2.3.As mentioned before some of the markers disappear for small amounts of time, in those cases, they are simply excluded from the computation of the errors (no interpolation is performed). 74 experimental validation repetition i ¯ e¯ s Polyester 0.67 cm 1.09 cm Wool 0.51 cm 0.97 cm Denim 0.50 cm 1.00 cm Stiff-cotton 0.38 cm 0.74 cm A2 0.71 cm 1.08 cm A3 0.33 cm 0.83 cm Shake 0.57 cm 1.01 cm Twist 0.47 cm 0.89 cm Slow 0.40 cm 0.80 cm Fast 0.63 cm 1.11 cm Global 0.51 cm 0.96 cm Table 4.5: Mean absolute error ¯ e and the mean standard deviation ¯ s with the optimal value of the parameters averaged over: fabric’s material, size (A2or A3), type of movement and speed for the first repetition of the recordings. Figure 4.15: Surface plot of the error function ¯ e(δ,α) for the fast shaking motion of A2wool (left) and a close-up near the detected minimum (right). Notice the presence of noise in the close-up. To find the optimal value of the parameters, we minimize the average on time of the absolute error (4.3) by performing a sweep search on the space (δ,α)∈[0, ρ]×[0, 4ρ] . The lower and upper limits are selected based on physical (they have to be positive) and empirical (if they are too large they drag the cloths too much, as if they were underwater) considerations. We found this optimization method to be faster and more robust that the more sophisticated minimization algorithm used in Section 4.1. This is likely due to the fact that the function we are trying to minimize is not completely smooth (see Figure 4.15) and hence has many local minima. The results for all 64 experiments can be checked in Appendix A. In Table 4.5we can see the mean absolute error and the mean standard deviation averaged over: material, size, 4.2 aerodynamics study 75 type of movement and speed for the first repetition of recordings. In Figure 4.16,4.17 and the video at https://youtu.be/UWda9cw0uI4 , we can see a visual comparison between the recording and its simulation with the optimal value of the parameters. Figure 4.16: Three frames comparing the recorded fast twisting of A2 polyester (left) with its inextensible simulation (right). The error at the three depicted frames from left to right is 1.82,1.73 and 1.36 cm respectively; being the average error of the whole simulation 1.15 cm. From the table, we can deduce that for the inextensible model, the most challenging material to modelize is polyester. This is somehow to be expected because of its silk-like properties. On the other hand, the errors are larger for the bigger textiles, this is again reasonable since they have double the area. The fast motions have larger errors, this is likely due to the fact that in that case aerodynamics are harder to model. Finally, we can see that both the shaking and twisting motions have comparable errors. Figure 4.17: Three frames comparing the recorded fast shaking of A2denim (left) with its inextensible simulation (right). The error at the three depicted frames from left to right is 0.63,0.92 and 1.001 cm respectively; being the average error of the whole simulation 0.84 cm. 4.2.3A priori forecast of αand δ As we can see in Tables A.1and A.2in Appendix A, the values of α and δ that minimize the absolute error ¯ e vary substantially with respect to material, size, speed, etc. We would like to find an a priori formula that we can use to forecast the value of these parameters without 76 experimental validation minimizing ¯ e(α,δ) . Apart from intrinsic properties of the textile (the density ρ and its size), this formula will also depend on the cloth’s speed. We look then for formulas of the form: δ=δ0+δ1S+δ2V+δ3ρ,α=α0+α1S+α2V+α3ρ, (4.6) where S is a measure of size and V a measure of velocity. We have found that for our purposes using as S a normalized area of the cloth ( 1 for size A3and 2 for A2) and as [V] = m2·s−2 the average (in time and over all the nodes) of 50% of the highest squared velocities gives the best results. Now, in order to find the optimal values of {δk,αk} we minimize the function given by: R(δ0, . . . , δ3,α0, . . . , α3) = 1 32 32 ∑ m=1 ¯ em(δ,α), (4.7) where ¯ em corresponds to the mean on time of the absolute error (4.3) of the m -recording of either repetition I or II obtained using the values α,δ given by Equation (4.6). To minimize this function we employ again a derivative-free algorithm (the Nelder-Mead Simplex Method). We denote by R∗ Function (4.7) evaluated at the optimal values of {δk,αk} . For comparison let E∗ be the mean of the optimal errors found in the previous section (which will be by definition smaller), i.e. the mean of the errors displayed in Table A.1or A.2. repetition i repetition ii δ0-0.0223 -0.0359 δ1-0.0178 -0.0117 δ20.0714 0.0780 δ30.7664 0.7890 α00.2082 0.2155 α1-0.1481 -0.1711 α21.1804 1.4410 α31.7440 1.9387 R∗0.578 cm 0.646 cm E∗0.516 cm 0.560 cm Table 4.6: Optimal values of the parameters δk,αkfor both repetition I (with hanger) and II (with bare hands) obtained by minimizing the function R. In Table 4.6we can see the optimal values of the parameters {δk,αk} along with the mean absolute errors averaged over all recordings for 4.3 validation of collisions 77 both repetition I (with hanger) and II (with bare hands). Notice that the error R∗ is comparable to E∗ and hence the fitting is quite accurate. Moreover, both sets of parameters are estimated independently for the two repetitions and have very similar values and the same signs, which shows that they are significant and have a consistent meaning. In particular, this justifies the introduction of the novel aerodynamic parameter δdone in Section 2.5. 4.2.4Discussion of the results In this second round of experiments, we have performed a more exhaustive set of recordings than before. We already knew from Section 4.1that our model was capable of reproducing faithfully the shaking motion of A3textiles, but we have enlarged the set of recordings by adding a new twisting movement and a larger cloth size (DIN A2). This has resulted in two sets of 32 recordings (each cloth being recorded two times on different days: with a special hanger or bare hands) each lasting approximately 15 seconds. As before we have estimated the optimal values of the physical parameters and the model has achieved very low mean errors (less than 1cm) and standard deviations (see Tables 4.5,A.1and A.2), even for the A2textiles and fast motions. Finally, we have found a predictive formula in order to obtain a priori estimates for the values of the parameters α and δ of the model. This formula depends on the density of the textile, its size and more importantly its speed. The formula was found to produce parameter’s values which in turn still give rise to very low absolute errors. In the next section, we will further make use of this a priori formula and test its accuracy. 4.3 validation of collisions For this final set of experiments, we will use again the following materials: polyester, wool, denim and stiff-cotton. We employ exactly the same recording setup described in the previous section (the motion capture system, Section 4.2). This time only the A2size is used (as before with 20 reflective markers placed equidistantly). We denote the sequence of positions of the recorded fabric’s nodes by {ϕ0,ϕ1, . . . , ϕm} and the simulated sequence by {φ0,φ1, . . . ,φm} . For the simulations, we utilize a refined 7×9 mesh. As error metrics, we use again the absolute error (4.3) and the spatial standard deviation (4.5). In order to validate the realism of our collision model, we must fit this time three parameters: as before α (damping) and δ (virtual mass); and for the first time µ (friction coefficient). In order to obtain their optimal value, as usual, we minimize the mean on time of the absolute error: 78 experimental validation ∑ i ei(δ,α,µ) = ∑ iq||φi(δ,α,µ)−ϕi||2 M. (4.8) The experiments are performed by a human (without any hanger) and consist of two scenarios: Figure 4.18: Putting a tablecloth motion sequence (right to left): the cloth starts suspended and is afterwards laid dynamically (only partially) onto the table. 4.3.1Tablecloth scenario The textile starts suspended at about 10 cm of height and is afterwards laid dynamically (only partially, so that half of the cloth is still suspended) onto the table (see Figure 4.18). Each motion lasts approximately 4seconds (with a frame every dt =0.01 seconds) and is performed with two different surfaces as the table, one with low friction (a raw polished table) and one with high friction (a table with a tablecloth). The goal here is to estimate the friction coefficient µ (see Equation (3.6)) for the two different surfaces and to study the sensitivity of the model with respect to friction. material µlow ¯ elow ¯ slow µhigh ¯ ehigh ¯ shigh Polyester 0 0.95 cm 1.20 cm 1 0.84 cm 1.03 cm Wool 0 0.58 cm 0.73 cm 2 0.52 cm 0.75 cm Stiff-cotton 0 0.60 cm 0.86 cm 2 0.58 cm 0.77 cm Denim 0 0.77 cm 1.11 cm 1.6 0.61 cm 0.80 cm Table 4.7: Optimal values of the friction coefficients along with the mean absolute error and spatial standard deviation for the low and high friction scenario In Table 4.7we can see the optimal values of the friction coefficients along with their optimal errors and deviations for the low and high 4.3 validation of collisions 79 friction scenarios. The optimal friction coefficients for the low friction case (raw polished table) were all smaller than 10−3 and that is why they were rounded up to zero on the table. For a visual comparison of the results see Figure 4.19 and the video at https://youtu.be/sWJcx fTwKHE. Figure 4.19: Three frames comparing the recorded tablecloth low friction scenario of A2wool (left) with its inextensible simulation (right). The error at the three depicted frames from right to left is 0.80, 1.17 and 0.76 cm respectively; being the average error of the whole simulation 0.58 cm. In order to understand how friction influences the dynamics of the textiles we perform a sensitivity analysis for the high friction case, i.e. we vary the value of µ (keeping all the other parameters fixed), and compute the mean of the absolute error (4.3). The results can be seen in the heat-map depicted in Figure 4.20. Notice that in general the model is quite stable with respect to the optimal friction value. Figure 4.20: Sensitivity analysis for high friction case, i.e. we vary the value of µ and compute the absolute error (4.3) for the four A2fabrics. In red we encircle the error found with the optimal parameter of µ. 80 experimental validation 4.3.2Hitting scenario In this final scenario, the fabrics are held suspended in the air (with the long sides perpendicular to the floor) and hit repeatedly with a long stick. The hits are aimed at various locations of the cloth with varied strengths and speeds (see Figure 4.21). In order to simulate the hits, the stick is subdivided into small edges and we employ a procedure similar to the one used for self-collisions of the cloth in the case of an edge-edge collision (see Section 3.3). Remark 4.3.1.The stick could be represented as a moving cylinder of small radius, but then we would be forced to use a very fine mesh to simulate the hits. Figure 4.21: Long-stick hits sequence (left to right): the cloth is held by its two upper corners and then is hit repeatedly with a long stick. The hits are aimed at different locations with varied intensities. Let us denote by {a1(t), . . . , am(t)} the endpoints of the edges of the stick. Then, for every iteration j of the iterative process (3.4), as we do with self-collisions (see Section 3.4.1), we must check if a collision occurred during the motion between the stick-edges {a1(tn), . . . , am(tn)} →dt {a1(tn+1), . . . , am(tn+1)} and the edges of our triangulated cloth φn→dt φj . This means that for every detected collision, in the next iteration j+1 of the sequence of quadratic problems (3.4) we must add a constraint of the form: H(φj+1) = ⟨πα(a1,a2)−πβ(x1,x2),ν⟩ ≥ 0, where a1,a2 are the two endpoints of the corresponding edge of the stick, x1,x2 are likewise the two endpoints of the edge’s cloth, πα′(a1,a2) = (1−α′)a1+α′a2 and πβ′(x1,x2)=(1−β′)x3+β′x4 are the closest points between the two segments and ν is the normal vector to both edges. The values ν,α′,β′ are constant in time, and are computed in the case of the cloth with the positions of the segments defined by φj and for the stick at time tn+1 . The normal vector ν is oriented such that H(φn)≥0. 4.3 validation of collisions 81 On the other hand, the real long stick has a length of 75 cm and a diameter of 1.5cm (see Figure 4.22). Two (rather large, with a diameter of 1.5cm) markers are put at both ends of the stick to record its trajectory. Figure 4.22: Long-stick (with a length of 75 cm and a diameter of 1.5cm) used to hit the textiles. Two markers are put at both ends of the stick to record its trajectory. Remark 4.3.2.We consider the stick’s thickness by imposing H(φ)≥ τ0 , where τ0=0.75cm is the radius of the stick. Moreover, this thickness is taken into account in the detection process (see Section 3.3.3). The stick is made of polished plastic and hence we consider friction between the cloth and the stick to be negligible (moreover, since the cloth is held firmly by the two upper corners the small amount of friction that could exist is always overcome by the stick). Each textile is hit four times with recordings varying between 12 and 18 seconds (as usual with a frame every dt =0.01 seconds). On top of fitting as usual the damping parameter α and the (virtual) gravitational mass δ , the goal of this scenario is to assess the realism of our collision algorithm when modeling the hits and put to use the predictive aerodynamic formula (Section 4.2.3). Remark 4.3.3.By the nature of this collision experiment, some movements of the textiles are very abrupt and therefore as mentioned before the markers disappear some of the time. This problem is also present in the recording of the trajectory of the stick. This is more problematic, since we need its position at every time to be able to simulate the collisions. We have interpolated the missing positions of the stick linearly. In this experiment, we also study the performance of the activeset collision algorithm described in Chapter 3. We compare it with a standard interior-point algorithm implemented to solve quadratic problems (see [86]). As before, all comparisons are performed using an Intel Core i7-8700K with 12 cores of 3.70 GHz. Since the four recordings have different durations, we compute the quotient 88 surface reconstruction via morse theory of sampling density or regularity, or of surface topology. It can be applied to surfaces in any ambient dimension. We use the gradient flows of [39,126] as the starting point, but then detect saddle points and their Morse cells differently, proposing a new procedure based on studying the level sections of these flows. 5.2 morse theory for manifolds Let M be a smooth compact manifold without boundary. A map f: M→R is Morse if it is C2 , has only finitely many critical points, and at all of these the Hessian H(f) is nondegenerate. Classical Morse theory (see [52]) shows that a generic Morse function f induces, through its gradient flow, two decompositions of the manifold M: 1. As a CW complex (see [83]): Each critical point of f , together with its unstable manifold for the vector field −∇ f , forms a cell which is topologically a ball, whose boundary attaches to lower-dimensional cells (see Figure 5.1). A global piecewise parametrization of M is achieved, and a Morse-Smale complex, with the critical points of f as a basis, giving the singular homology of M. Figure 5.1: Critical points of the Morse-Smale function f(x,y,z) = z on an example surface. 2. As level sets: M is foliated by the level sets f−1(c) . For regular values c these level sets are submanifolds of M with codimension 1, with f−1(c1)∼ =f−1(c2) if no critical value of f lies between c1 and c2 . The transformation of the level set when c crosses a critical value of fis a surgery (see [52]). The success of Morse theory comes from the fact that Morse functions, and the Morse-Smale transversality conditions required for the above analysis, are generic among C2 maps from M to R . For instance, the height function in a random direction in RN has probability 1 of being a Morse-Smale function. In the following subsection, we present a summary of Morse theory (as level sets) for surfaces without boundary with views towards explaining the case with boundary next. 5.2 morse theory for manifolds 89 5.2.1Morse theory for surfaces without boundary Let S⊂RN be a smooth compact surface without boundary. As before, a map f:S→R is Morse if it is C2 , has only finitely many critical points (i.e. points p where dpf=∇f(p) = 0 ), and at all of these its Hessian has rank 2. Definition 7(Morse data).For each critical point p∈S , the Morse data are the pair of sets (A(p),B(p))where A(p):=B(p)∩f−1([f(p)−ϵ,f(p) + ϵ]), with B(p)⊂RN a closed ball around p of sufficiently small radius, the value of ϵ>0 is such that there are not more critical points of f in f−1([f(p)−ϵ,f(p) + ϵ])and B(p):=B(p)∩f−1(f(p)−ϵ). Notice that B(p)⊆∂A(p). Let us now denote f≤c=f−1((−∞,c]),fc=f−1(c). Then it is well known (see [52]) that as c∈R increases two things can happen: a If c1<c2 and there are no critical values of f between them, then f≤c1 and f≤c2 have the same topology (they are actually diffeomorphic). b If there is a critical point p∈S such that c1<f(p)<c2 then f≤c2 is obtained, up to diffeomorphism, from f≤c1 by attaching the cell A(p) along B(p) , i.e. f≤c2≃(f≤c1∪A(p))/∼ where the equivalence relation is given by identifying B(p)⊆f≤c1 with points of ∂A(p)⊇B(p). Now, since H(f) has rank 2at p , we can only have three types of critical points (see Figure 5.2) according to their index (number of negative eigenvalues): 1. Minima: A(p) is homeomorphic to a disk and B(p) = ∅ . In this case, there is no surgery, the cell A(p) just appears. This cell retracts to a point, the local minimum, which will be a 0 -cell in the Morse-Smale complex. 2. Saddles: A(p) is a quadrilateral (homeomorphic to a disk) and B(p) consists of two segments. Two opposite sides of ∂A(p) are identified with B(p) (see Figure 5.2). The attachment of the cell A(p) along B(p) is homotopy-equivalent to the attachment of a 1 -cell, namely the medial axis of A(p) to the middle points of B(p). 90 surface reconstruction via morse theory 3. Maxima: A(p) = D is homeomorhic to a disk and B(p) = ∂D is its boundary. The attachment map identifies ∂A(p) with B(p) . This surgery adds a 2-cell to the complex. Figure 5.2: Three types of critical points (from top to bottom: minima, saddles and maxima) and their Morse data for surfaces without boundary. 5.2.2Morse theory for surfaces with boundary Morse theory also extends to manifolds with boundary [46]. Since this theory is less well known, in the following we give more details on how the boundaries affect the cellular decomposition and the transition between level sets. Definition 8.We say that a C2 map f:M→R defined in a manifold Mwith boundary ∂Mis Morse if 1. it is Morse in the interior of M, 2. its restriction to ∂M,g:=f|∂Mis also a Morse function, 3. if p∈∂Mis a critical point of gthen ker(dpf) = {0}. As in the previous section, we will focus on the case M=S is a surface. Notice that now we can have critical points of g in the boundary curves ∂Sthat are not critical points of fin the whole of S. Moreover, the third condition of the previous definition ensures that critical points located at ∂S are not saddle points of f in S . Notice that minima (resp. maxima) points p of f located at ∂S will satisfy that 5.2 morse theory for manifolds 91 dpf=0, and thus, they will not be strictly speaking critical points of f(but they will be of g). Hence, we will adopt a Whitney stratification dividing the surface into two strata: the first E1=∂S being the boundary and the second E2=int(S) the interior. Each stratum will have as before Morse data, but this time apart from the sets (A(p),B(p)) defined before, which we will call tangential data, we will have also normal data. Definition 9(Tangential and normal data).Let p∈S be a critical point for some stratum Eiof the Morse function f|Ei. Then 1. The tangential Morse data are the pair of sets (AT,BT) defined as in Definition 7for f|Ei. 2. The normal Morse data are the sets (AN,BN)given by AN(p):=N(p)∩f−1([f(p)−ϵ,f(p) + ϵ]), BN(p):=N(p)∩f−1(f(p)−ϵ). where N(p) = S∩D(p) is called the normal slice at p and D(p)⊂ RN is a sufficiently small closed disk around p of dimension N−dim(Ei) transversal to Ei at p and the value of ϵ>0 is such that there are no more critical points in f−1[f(p)−ϵ,f(p) + ϵ]. Notice that by construction Ei∩D(p) = {p} , and as before BN,T⊆ ∂AN,T. Remark 5.2.1.When p∈int(S) and N=3 then D(p) is homeomorphic to a piece of curve normal to S at p and hence N(p) = {p} . Therefore (AN,BN) = (p,∅) . It is not hard to see that this is also the case when N>3. We are now ready to state the main theorem of stratified Morse theory [46] (SMT theorem, pages 6-8). Theorem 1(Goresky-MacPherson).As c∈R increases two things can happen: a If between c1<c2 there are no critical points of f then f≤c1 and f≤c2are diffeomorphic. b If there is a critical point p∈S such that c1<f(p)<c2 then f≤c2 is obtained from f≤c1 by performing a surgery around p with Morse data (A,B)diffeomorphic to the topological product (AN,BN)×(AT,BT) = (AN×AT,AN×BT∪BN×AT), i.e. f≤c2≃(f≤c1∪A)/∼ where the equivalence relation is given by identifying B⊆f≤c1with points of ∂A⊇B. 92 surface reconstruction via morse theory Remark 5.2.2.We already established that at interior critical points of S , we have (AN,BN) = (p,∅) . Therefore in that case the above product is trivial and the attachment maps are the same ones described earlier. The new cases occur when p is a critical point of f|∂S lying on ∂S . Since ∂S is a one-dimensional curve, p∈∂S can only be a maximum or a minimum. Nevertheless, depending on whether p is also a local minimum (resp. maximum) of f or just of g=f|∂S , we will have four different cases (see Figure 5.3). Recall that for minima (and maxima) points p of f located at ∂S we have that dpf=0 . Nevertheless, in order not to complicate the discussion semantically we will still call them critical points since they satisfy dpg=0. Remark 5.2.3.We will denote by (■,⊔,| |, _) a quadrangular 2 -cell, all its sides minus the top one, two lateral sides and the bottom side, respectively. Figure 5.3: Four types of critical points located at the boundary and their tangential and normal Morse data. The four new different critical points that we can have are: Maxima of f|∂S : In this case, AT is a closed concave piece of curve (which we denote by ∩ ) and BT two points, i.e. (AT,BT) = (∩, . .) . We now have two sub-cases: 1.p is also a local maximum of f . Then AN is a closed interval and BT one point, i.e. (AN,BN) = (|, . ) . Therefore the topological product is (A,B) = (∩, . .)×(|, . )=(■,⊔), 5.2 morse theory for manifolds 93 and the three of the sides of ∂A where p is not present are attached to B . During this surgery, a 1 -cell (containing p ) is attached to the boundary and a 2 -cell to the interior (closing a void in the process). 2.p is not a local maximum of f (only of f|∂S ). Then AN is again a closed interval but BT is empty, i.e. (AN,BN) = (|,∅). Then the Morse data is (A,B) = (∩, . .)×(|,∅)=(■,| | ), and the two opposite sides of ∂A where p is not present, attach to B . This surgery is equivalent to attaching a 1 -cell to the boundary (but no 2-cell is attached to the interior). Minima of f|∂S : In this case, AT is a closed convex piece of curve (which we denote by ∪ ) and BT is empty, i.e. (AT,BT) = (∪,∅) . We again have two sub-cases: 1.p is also a local minimum of f . Then AN is a closed interval and BT is empty, i.e. (AN,BN) = (|,∅) . Therefore the topological product is (A,B) = (∪,∅)×(|,∅)=(■,∅). In this case there is no surgery, the cell A just appears. This cell retracts to a point, which will be a 0 -cell in the cell complex. 2.p is not a local minimum of f (only of f|∂S ). Then AN is a closed interval and BT is a point, i.e. (AN,BN) = (|, . ) . The Morse data is (A,B) = (∪,∅)×(|, . )=(■, _). and the opposite side of ∂A where p is, attaches to B . This surgery is equivalent to attaching a 0 -cell to the boundary and a 1-cell to the interior. We now make a summary of the effect of each critical point on the cell complex once we have made the appropriate deformation retracts. 1. Interior maximum: attach a 2 -cell to the 1 -cells or to the boundary curves. 2. Interior saddle: attach a 1 -cell to the 0 -cells or to a point of the boundary. 3. Interior minimum: add a 0-cell (a point) to the skeleton. 4. Local maximum of S located in ∂S : attach a 2 -cell to the 1 -cells and attach a 1-cell to a point of the boundary. 94 surface reconstruction via morse theory 5. Boundary maximum (not of S ): attach a 1 -cell to a point of the boundary. 6. Local minimum of S located in ∂S : add a 0 -cell to the skeleton (this point will be on the boundary). 7. Boundary minimum (not of S ): add a 0 -cell (boundary point) to the skeleton and attach a 1 -cell to a local minimum ( 0 -cell) or to a point of the boundary. In order to construct the complex, we first add all the 0 -cells (all interior and boundary minima and possibly some boundary points), then the 1 -cells (the boundary curves, the 1 -cells corresponding to each saddle point and the 1 -cells joining boundary minima to local minima or boundary points) and finally we add the 2 -cells corresponding to each local maximum. 5.3 extension to point-clouds In this section, we outline the main ideas of the reconstruction algorithm for point-clouds. First, we explain the case without boundary and then with boundary; in the next section, we will give full details for both cases. Let X⊂RN be a point cloud sampling a compact surface S , possibly with boundary. Determine the neighbors of each point p in the sample, e.g. proceeding as outlined in Section 5.4. Choose a unit vector ν∈RN such that the height function f(p) = p·ν has different values in all points of X , and define the positive gradient (or upwards), resp. negative gradient (or downwards) flows of f by sending every point p in the cloud to its neighbor that maximizes the slope of growth of f , resp. makes f decrease with the most negative slope. Points where f cannot grow, resp. decrease, are local maxima, resp. minima of f in X . The delicate critical points to compute are saddles. Interpret the downwards flow of f on X as an embedded graph whose vertices are the points in the cloud and the edges are given by connecting each point in the cloud to its downwards neighbor. The intersections of this graph with the hyperplane x·ν=c can be considered point samples for the level set f−1(c) on the original surface S , with some noise added by the linear interpolation. This level set consists of point samples of curves, either closed or with edges in the boundary of S . These level curves can be reconstructed, e.g. as explained in Section 5.4, identifying the different components. Perform these level set intersections at n different levels ci=c0+i·h ranging from c0=min f(X) to cn=max f(X) . The number of level sections n must be selected, the idea is that all surface features (e.g. saddle points) whose range in height is h or greater will be detected. 5.3 extension to point-clouds 95 In Section 5.4we will explain how changes in the topology of the level set correspond to variations in the number or type of curves. For example, going over a saddle point of f can be tracked by detecting two pairs of neighbors in different connected components of a level set which end up in the same connected component of the next level set after applying the downwards flow (see Figure 5.4). Figure 5.4: Change in level set when crossing a critical value in a surface (left) and point cloud (right): note the change in neighbors among the 4marked points after the flow. Following the 4points in these 2pairs in the downwards flow, and their pairing according to closeness, the level at which the pairing changes marks the position of the saddle point (see Figure 5.4). The unstable variety of this saddle point for the flow of −∇ f (i.e. the 1 -cell joining the saddle point to local minima) is approximated by taking the two pairs of neighboring points at the level of f immediately bellow the saddle point ( {B′,D′} and {A′,C′} in the figure), and averaging the downward orbits of each pair (which end up in a minimum, but not necessarily the same). These computations have a margin of error O(d) , where d is the local variation of height among neighboring cloud points. Once the saddle points of the height function and their (un)stable varieties have been found, the piecewise parametrization for the entire surface without boundary follows: 0 -cells are the local minima, 1 -cells have been parametrized at the saddle point detection, and each 2 -cell can be identified from the tree formed in itself by the upwards flow to its unique maximum. Finally, the boundary relations given by the downwards flow on the cells give us the Morse-Smale complex and singular homology of the surface S. The case of surfaces with boundaries poses several complications that we will discuss in detail in the next section, but the main aspect to take into account with respect to 1 -cells is that when we have boundary minima (i.e. a critical point of f|∂S that is a minimum when restricted to the boundary but not in the whole point-cloud) we must add to the cellular decomposition of S a 1 -cell that joins the boundary minimum with a local minimum or to a boundary point of the surface. This is 96 surface reconstruction via morse theory simply obtained by applying the downwards flow repeatedly starting at the boundary minimum. 5.4 practical implementation In this section, we give detailed algorithms for all the ideas presented in the previous section. 5.4.1Neighbors identification The first step is the identification of a set of neighbors of each point v in the cloud X⊂RN. There are two classical approaches: 1. k-nearest neighbors (KNN): given a value for k and a point v , the k nearest points {v1, . . . vk} with respect to the Euclidean distance are declared as its neighbors. This is quite efficient to compute but runs into problems when the point-cloud has irregular densities and k is not big enough, e.g. when all the closest points to v are clustered at one side of it and do not enclose the vertex (see Figure 5.5, left). 2. Voronoi-Delaunay neighbors: we perform Voronoi’s cellular decomposition of the point-cloud and then declare as neighbors of v the points vi with neighboring cells (i.e. are connected to it by an edge in the Delaunay triangulation). This has the virtue of enclosing the vertex v even with irregular densities, but it can be expensive to compute and produces neighbors which are too apart from each other (see Figure 5.5, right). Figure 5.5: Typical problems associated to k-nearest neighbors (left) and Voronoi neighbors (right). On the left, the majority of the closest points to v are clustered at one side of it. On the right, vertices that are too far apart from each other have a neighboring cell. Therefore we merge these two criteria, and declare two points as neighbors when (i) each point is among the k -nearest neighbors of the other and (ii) their Voronoi cells in the decomposition of the ambient space RN induced by X are adjoining. In order to be efficient we first choose a k≈12 to compute the k -nearest points and then 5.4 practical implementation 97 only keep as neighbors the vertices that are connected to v in the Delaunay triangulation of these few points. Finally, the relationship of neighborhood is made symmetric by reciprocating neighboring relationships where needed. The neighbors of v will be denoted by Neigh(v) . This produces a (locally non-planar) graph, which gives an idea of the local structure of X , but which will be in general very complicated. Remark 5.4.1(Pruning of the point-cloud).If the original cloud X has a very irregular density, it can be wise to discard some points. In order to do this, several heuristics can be employed, one of them is to estimate the density at each point (e.g. the reciprocal of the logarithm of the distance to the closest point) and then discard points whose density is less than the mean plus a fixed multiple of a standard deviation. Special care must be taken to avoid removing complete clusters of points (at least one should be kept). 5.4.2Tangent space estimation This task is performed through Principal Component Analysis: if the point v and all its neighbors Neigh(v) = {v1, . . . vk} were co-planar we would have that for every i : ⟨ vvi,nj⟩=0 , where nj are all the normal vectors to the surface (recall that in general we are in RN ). Since in general this will not be the case, we find the nj ’s by minimizing the function ∑j∑k i=1⟨ vvi,nj⟩2 . This is equivalent to finding the regression plane in the least squares sense, and it can be done efficiently by means of a singular value decomposition of the matrix with vectors  vvi as rows. 5.4.3Boundary recognition Once we have an estimation of the tangent spaces, in principle a boundary point of the surface can be easily identified because after orthogonally projecting it and its neighbors on its tangent plane, they cluster in a semi-space (see the first panel of Figure 5.6). Nevertheless, this method is difficult to implement robustly (e.g. on points of high curvature of the boundary curves, see the second panel of Figure 5.6). In order to obtain a robust detection, the idea will be to declare points as lying on the boundary only when the projections do not enclose the point. In points with high curvature where the tangent plane may not be perfectly estimated the previous method can give false positives. In order to overcome this difficulty, we will also project the point and all its neighbors in the tangent planes estimated for the neighbors. We will build a graph for every projection and only declare v as boundary point when none of the graphs enclose v. 104 surface reconstruction via morse theory Figure 5.13: A sampled dumbbell: the black line is the direction of the height function; local maxima, resp. minima, are painted red, resp. black; saddle points are painted blue; their 1–cells are outlined in blue. There are two 2 -cells: one in dark red (left) and one in light blue (right). 1. Boundaries: when S has non-empty boundary, the boundary curves are also part of the 1 -skeleton. We already parametrized those. Each boundary maximum gives rise to a 1 -cell that attaches to a 0-cell (point) located at the boundaries. 2. Saddle points: there are one-dimensional curves that go from one local minimum m1 or boundary point to another (not necessarily distinct) local minimum m2 or boundary point passing through the saddle point. 3. Boundary minima: there are 1 -cells in the complex not associated with a saddle point or to boundary curves. In this last case, they always connect a boundary minimum to a local minimum or boundary point of the cloud. In order to compute this 1 -cell we simply flow down every boundary minimum using Down(·). Remark 5.4.5.When 1 -cells introduced by saddle points or boundary minima end at points of the boundary ∂S not previously labeled as 0 -cells, we introduce the point where they meet to the 0 -skeleton and divide the 1-cell into two. Computation of saddles and their 1-cells We now explain how to compute saddle points and afterwards their associated 1 -cells. Recall that we know when we are in the presence of a saddle: when going from Γ(ci+1) to Γ(ci) a connected component appears, disappears or there is a change of bordism pairing of boundary points. In any case, all points of Γ(ci+1) have images in Γ(ci) using the flow Down(·) . Therefore, there exist four points v1,v2,v3,v4 in Γ(ci+1) that can be paired by proximity (i.e. are neighbors in the level 5.4 practical implementation 105 set curves) as {v1,v2} and {v3,v4} whose images by Down(vi) = ˜ vi are now paired differently as {˜ v1,˜ v3} and {˜ v2,˜ v4} . The saddle point is approximated by the average: s=1 8 4 ∑ i=1 vi+ 4 ∑ i=1 ˜ vi!. One branch of the 1 -cell is obtained by taking the average of the two orbits generated by Down(·) starting from {˜ v1,˜ v3} . The other branch is approximated in the same manner but starting from and {˜ v2,˜ v4} . This averaging operation generates new points not previously in the cloud. All these new points are added to the cloud, declaring as neighbors of s the 8points used for its computation, and for each new point of the branches its two neighbors in the 1 -cell plus the two points that were used for its computation. The flows Down(·) and Up(·) are simply defined at those points by following the newly created trajectories down or up (except the saddle points which are fixed points of both flows). This ensures that this newly defined curve (the 1-cell) is invariant by both flows. 2-cells There is one 2 -cell for each local maximum of f (recall that boundary maxima only generate 1 -cells, i.e. the boundary curves), we already found them (they are the fixed points of Up(·) ), but we will now explain how to deduce which points of the cloud correspond to each 2 -cell (or each maximum). The idea is simply to flow up every point v∈X by Up(·) and see at which maximum it ends up. In symbols, this means that we consider the limit of the sequence vn=Up(vn−1) where v0=v when n→+∞ (this limit exists because we have a finite number of points and maxima are fixed by Up(·) ). It can happen that the sequence {vn}n converges to a saddle s . This can only happen for points of the 1 -cells, but must be corrected for the rest: when a point v ends at a saddle s but it is not part of the 1 -cell, we look at the neighbors Neigh(v) and see in which maximum they finished. The most common maximum among the neighbors is selected as the corrected destination of v . It may be necessary to iterate this procedure a finite number of times. 5.4.8Attachment maps of the Morse cells Now that we have the skeleton of S , i.e. the 0, 1, 2 -cells, we encounter a delicate problem: figuring out the attachment maps between the cells. We already know how the 1 -cells attach to the 0 -cells. We now discuss how 2 -cells attach to 1 -cells. We must figure out which 1 -cells are the boundary of the 2 -cells, in which order, orientation and how many times they appear: once (when they are part of the boundary of 106 surface reconstruction via morse theory S or when they attach to a different 2 -cell) or twice (when they will be, formally, two different sides of the 2 -cell that are identified, see Figure 5.13). We first discuss the case without boundary. Smooth case We first describe the process for the smooth case (without boundary) and then explain how to adapt it to point-clouds. Notice that the 2 -cell corresponding to each maximum ˆ m , let us call it Dˆ m , is the set of points that converge to ˆ m under the flow +∇f (in Dynamical System terminology, its stable manifold). The main idea to deduce how the boundary of Dˆ m attaches to the 1 -cells is to choose a simple closed curve around each maximum and flow it down by −∇ f until it reaches the 1 -cells. To achieve this properly, we must perturb f so that in a neighborhood of the 1 -cells the gradient flow is transversal (i.e. not parallel) to the 1 -cells. This is needed since otherwise the points of Dˆ m would converge to minima under the downwards flow −∇ f . In order to obtain this perturbed ˜ f we consider an arbitrary small tubular neighborhood around the 1 -cells and a normal vector field to the 1 -cells. These two flows are joined smoothly in the tubular neighborhood with the aid of a partition of unity (i.e. bump functions, see [52]). Discrete case In the discrete case, we choose a simple (i.e. without self-intersections) closed curve of neighbors around each maximum (a crown) and flow it by Down(·) until it reaches the 1 -cells. We do this in a way such that when a point of the curve has neighbors that are in the 1 -cell, we stop following the flow Down(·) and match the point in question to the point in the 1 -cell with results in the most negative downwards slope. This method would work perfectly if simple closed curves were preserved by Down(·) ; since this is not the case, we must repair the curve every time we flow it down. There are three main situations: 1. Splitting of points: it may happen that two points w1 and w2 were neighbors but Down(w1) and Down(w2) are not. In that case, we define the sub-graph of neighbors given by points {v∈X:f(v)≤max [f(Down(w1)),f(Down(w2))]}, and compute the shortest path joining Down(w1) to Down(w2) (taking as the distance of a path, that given by the edges of the graph). 2. Collapse of points: it may happen that two points w1 and w2 have the same image Down(w1) = Down(w2) . In that case, we simply delete one of the repeated points from the curve. 5.4 practical implementation 107 3. Creation of spikes: it may happen that one of two neighboring points w1 and w2 has image Down(w2) = w1 . This could happen to more than one point at a time. In that case, we simply delete the spiky points (i.e. Down(w2) = w1) in the new curve. The final problem we face is that the flowing crown may arrive at a stage where it is stationary by the downward flow but some points of it have not reached the 1 -cells. In that case, we pair each unmatched point of the crown with the closest (as measured by the edges of the neighboring relationship of the cloud) point of the 1-cells. We now have a matching between the initial crown and the 1 -cells (although not a bijection). But this is enough to deduce which 1 -cells are the boundary of Dˆ m (only the 1 -cells that are reached), in which order (this can be deduced since the crown is parametrized), the orientation (again thanks to the parametrization) and how many times they appear (once or twice, depending on the matching). Attachment maps for surfaces with boundaries When the surface S has a boundary ∂S the boundaries of the 2 -cells are no longer only the 1 -cells derived from the saddle points, but also possibly curves of ∂S . We distinguish several cases depending on the type of maximum we have: 1. Interior maxima (not belonging to ∂S ): as before we flow down a crown around the maximum until it reaches a point that is either a neighbor of a boundary point or of a 1-cell. In this way, we obtain a matching between the boundary of the 2 -cell and the 1-cells and boundary curves of S. 2. Boundary maxima (but not a local maximum of S): they do not intervene in the attachment maps of 2 -cells to 1 -cells (although they give rise to a 1-cell in the boundary). 3. Local maxima (located at ∂S ): we consider a semi-crown around the maximum ˆ m , meaning that we take an arc of the boundary centered at ˆ m and complete it with neighboring interior nodes so that we obtain a closed simple curve. Then we flow down this semi-crown in such a way that the arcs of the boundary stay in the boundary (we use the restricted flow ∂Down(·) ) and the interior nodes of the curve flow down normally by Down(·) . As before we flow down the semi-crown until each point of it reaches a node that is either a neighbor of a boundary point or of a 1-cell. For an example of this process in action, see the video at https: //youtu.be/8cgr54oRf6w. 108 surface reconstruction via morse theory 5.4.9Parametrization of the 2-cells Once we have each Morse cell and their attachment maps identified, the last problem we face is how to parametrize the 2 -cells, i.e. finding a flat domain D⊂R2 and a map ϕ:D→X , such that each ϕ(D)⊂X corresponds to one of the 2 -cells found before (see Figure 5.14). We will further require D to be a convex polygon and that ∂D is isometric to ϕ(∂D) . Once we have a parametrization of each 2 -cell, since they attach well (their boundaries are the 1 -cells), we have a full piece-wise parametrization of X. Remark 5.4.6.We may require ∂D to also have the same number of sides of ϕ(∂D) (i.e. the different 1 -cells), their length and their order (these are known in advance and will be determined by the attachment map of the 2 -cell to the 1 -cells). In that case, we denote the sides’ lengths by l1, . . . , ld. In order to find Dand ϕ, we follow two steps: Figure 5.14: Parametrization of a 2 -cell (left) using a rectangle D with isometric sides to the boundary of the 2-cell (right). 1. Obtain a convex polygon D⊂R2 whose sides are l1, . . . , ld : we find the polygon with maximal area and predetermined sides l1, . . . , ld inscribed in a circumference. To do so, we solve an optimization problem in order to find the circumference in question. By a result of Bramagupta [13], this is done by finding the radius that makes the polygon’s interior angles add up to 2π . Once we have D , we obtain a bijection between ∂D and the points of the corresponding 1-cells in X. Alternatively, we can take any convex polygon in the plane whose boundary ∂D is isometric to the boundary of the 2 -cell (e.g. a rectangle). 2. Obtain a correspondence of interior points of D mapping to the cloud points in the 2 -cell: for each interior point v in the 2 -cell we 5.4 practical implementation 109 want to find an interior point x=ϕ−1(v) in D . We assume that each point pi of D (including ∂D ) is a barycentric combination of its neighbors (as given by the neighboring relationship of the cloud, see [37]) with Tutte’s weights (all coefficients are equal to 1 |Neigh(v)| ). Then we have to solve a linear system with a unique solution (see [37]), which in turn gives us the coordinates of interior points of Dmapping to the cloud. Once we have these points, we can extend the parametrization to the whole interior of D by taking its Delaunay triangulation with the newly found vertices, and interpolating linearly for the images. Then we can re-mesh (or even quadrangulate) the polygon D if desired. This re-meshing allows us to obtain C∞ parametrizations of the surface via splines, to de-noise the point-sample if needed or to obtain a synthetic resampling of the point-cloud with better regularity properties. Remark 5.4.7.As already mentioned, the polygon does not need to have the same number of sides as the 2-cell, in fact, in order to apply some de-noising algorithms it is necessary to take the polygon D as a square or rectangle. Moreover, if the 2 -cell is too far from being flat (e.g. because of high curvature) the algorithm may produce interior points in D with non-uniform densities and thus it may be better to subdivide the cell into smaller and flatter pieces and afterwards parametrize each piece independently. On the other hand, using a polygon where the sides are mapped isometrically to each of the bounding 1 -cells, has the advantage of allowing the gluing of the resulting parametrizations of the 2 -cells into a global piece-wise defined and continuous parametrization of the entire surface. 5.4.10 Results In this section, we reconstruct three different surfaces (one without boundary –a torus– and two with it –a vest and a pair of pants–) which pose various challenges to the presented algorithm. The first two pointclouds are synthetic whereas the last one is a real scan of an actual textile. The three cases are challenging and interesting for different reasons: the torus is really slim and its embedding describes a (2,3)- toric knot which causes far away parts of the surface (as measured by geodesic distance) to be really near each other in euclidean space. The vest was obtained by cutting out parts of an ellipsoid in order to obtain a surface with the same topology as an open vest. Therefore, it has positive Gaussian curvature everywhere and very large boundary curves. Finally, the pants are a 3D-scan of a real pair of jeans and thus its point-cloud presents wrinkles, noise and an irregular density distribution of points. 110 surface reconstruction via morse theory Figure 5.15: A sampled knotted torus: the black line is the direction of the height function; local maxima, resp. minima, are painted red, resp. black; saddle points are painted blue; their 1 -cells are outlined in blue. On the bottom, we plot the level set curves highlighting when critical points appear. 5.4 practical implementation 111 Knotted torus Figure 5.15 shows our algorithm applied to a point-sample from a torus embedded in R3 along a (2,3)-toric knot. The algorithm correctly detects 2local maxima, 2local minima and 4saddle points for the height function depicted in the figure. Out of a point cloud of 30 000 points, a decomposition of the surface into 8Morse cells is found (two 0 -cells: the minima, four 1 -cells associated to the saddle points and two 2 -cells, one for each maximum). On the bottom of Figure 5.15 we display the level set curves Γ(c) = Gdown ∩Hc (see Section 5.4.6). Notice that for this point-cloud we have all possible local transformations of level-set curves for surfaces without boundary (see Figure 5.11). Ellipsoidal vest Figure 5.16 shows our algorithm applied to a sample of 36 000 points from a vest embedded in R3 . The cloud was obtained by cutting out parts of an ellipsoid in order to obtain a surface with the same topology as an open vest. After detecting and parametrizing the boundary successfully, the algorithm correctly detects 2local maxima, 2boundary maxima, 1local minimum and 3boundary minima for the height function depicted in the figure. All critical points are located at the boundary. Moreover, two new points (shown in purple) are added where the 1 -cells meet each other or the boundary curves (see Remark 5.4.5). Then, a decomposition of the surface into 16 Morse cells is found (five 0 -cells, nine 1 -cells and two 2 -cells). In order to deduce how the two 2 -cells attach and which 1 -cells are their boundary we apply the curve flow explained in Section 5.4.8. This process in action for one of the 2 -cells can be visualized at https://youtu.be/8cgr54oRf6w . Thus, we deduce how the cells attach with each other (e.g. the 1 -cell number 5appears on both 2 -cells and it is precisely one of the curves where they attach, see Figure 5.16). From this, we recover the entire topology of the vest. In Figure 5.17 we show a parametrization by a rectangle of one of the 2 -cells of the vest (the red one on the right in Figure 5.16). This is done as explained in Section 5.4.9: the bounding 1 -cells are mapped isometrically to a flat rectangle (notice that we consider the 1 -cells number 6and 6’ as different) and then interior points are obtained using the neighboring relationships of the cloud (taking care of removing neighbors at opposite sides of the 1 -cell number 6). Finally, the 2 -cell consisting of 27 000 points (shown in red in Figure 5.17) is down-sampled using the parametrization and interpolating linearly to a cloud of 900 points (shown in blue Figure 5.17). 112 surface reconstruction via morse theory Figure 5.16: A sampled vest: the black line is the direction of the height function; maxima are painted in red, minima in black; 1 -cells corresponding to boundary minima are outlined in blue and the boundary curves in black. The two purple points where the 1 -cells meet each other or the boundary curves are added to the decomposition. The numbers correspond to the different formal 1 -cells that, when identified (e.g. 7with 7’), reconstruct the entire surface from 2pieces homeomorphic to disks. 5.4 practical implementation 113 Figure 5.17: Parametrization by a rectangle of the rightmost 2 -cell of Figure 5.16. The bounding 1 -cells are mapped isometrically to a flat rectangle and then interior points are obtained using the neighboring relationships of the cloud. The 2 -cell consisting of 27 000 points (shown in red) is down-sampled interpolating linearly to a cloud of 900 points (shown in blue). 3D-scan of pants Figure 5.18 shows our algorithm applied to a 3D-scan of a pair of real jeans. The scan was made by a Artec Eva professional handheld 3D scanner while a person was wearing the garment. The 3D scan is processed with the proprietary software Artec Studio in order to obtain a point-cloud sample. In order to apply our algorithm and to reduce irrelevant details and noise, we down-sample the initial cloud of more than 500 000 points to 11 000 by using a box-grid filter, i.e. an axis-aligned bounding box is computed for the entire point-cloud and then divided into grid boxes. Points within each grid box are merged by averaging their locations. Still, the cloud has a lot of detail (e.g. wrinkles) that cause the proliferation of critical points and hence of Morse-cells. After computing for each point its neighbors as explained in Section 5.4.1, we apply a 120 developable surfaces Theorem 2(structure theorem for developable surfaces).A C3 developable surface S embedded in R3 has an open subset that is ruled, with unit normal vector constant along each line of the ruling but varying in a transverse direction. Every connected component of its complement is contained in a plane. This structure can be deduced from the Gauss map of the surface: Gaussian curvature 0 makes its rank 0 or 1 , the latter rank being reached on an open subset of the surface. The normal vector is locally constant in the rank 0 subset. The dichotomy in the rank of the Gauss map, and varied classical notations, motivate: Definition 12.A developable surface is torsal if the Gauss map has rank 1 on a dense open subset. Flat patches are connected subsets of a developable surface with nonempty interior where the Gauss map is constant, i.e. they are contained in a plane. The subdivision of a developable surface into torsal and flat patches is given by the boundary of the vanishing locus of the mean curvature (the trace of the Gauss map) and is not necessarily simple. Let us fix a planar domain R⊂R2 which is compact, contractible and has a piecewise C∞ boundary (e.g. a convex polygon). Define S to be the set of all C3 surfaces in R3 isometric to R . These surfaces are all developable, and S may be seen as the space of states of an inextensible (i.e. isometric for the inner distance) deformation of R in Euclidean space. The space of states S can also be defined as the set of C2 maps from R to R3 which are isometries with the image. As such, it is endowed with the compact-open topology derived from the Euclidean one in R and R3 . This topology furnishes valuable tips for the study of S : the set of surfaces containing flat patches has an empty interior because there exist arbitrarily small deformations making the normal vector nonconstant on an open set. Torsal surfaces are stably torsal if the mean curvature function intersects transversely the zero function. These surfaces form an open subset U ⊂ S , and suffice for our practical study of S. We can try to develop coordinates for the stably torsal state space U based on the classical structure theorem. First, let us recall how to identify developable surfaces among the ruled ones: Proposition 2(classical, see [19]).A ruled surface parametrized as ϕ(u,v) = γ(u) + v·w(u) , where γ is a regular parametrized curve and w a vector field over γ , is developable if and only if the 3vectors γ′(u),w(u),w′(u)are linearly dependent for all u. Given a regular C2 curve γ there is a way to obtain systematically such rulings over γresulting in regular torsal surfaces: Proposition 3.Let n be a unit normal C1 vector field over a regular, C2 curve γ with n′=0 . Then w=n×n′ defines a torsal surface in 6.1 the space of developable surfaces 121 a neighborhood of γ . Moreover, all regular, torsal rulings over γ are generated by such w, and only n,−ndefine the same torsal surface. Proof. n is normal to γ′ and w by their definitions, and w′=n×n′′ so at every u the vectors γ′,w,w′ are normal to n . Also, note that w=0 because n′=0 . If ˜ n is another unit normal vector field such that ˜ nט n′=µw for some function µ(u) then note that ˜ n has to be normal to both w and γ′ , hence a multiple of n . Finally, let us point out that if w is a nonvanishing tangent vector field over γ defining a torsal surface around it, then we can select a unit vector field n normal to γ′,w,w′ . The fact that n is normal to w and w′ imply that n′ is also normal to w, so n×n′is a multiple of w. To define coordinates in the space of stably torsal surfaces S isometric to a fixed bounded domain R , the pairs (γ,n) of Proposition 3run into a practical difficulty: the condition that n′=0 forces the Gauss map to have rank 1. If S is a stably torsal surface with mean curvature H of varying sign, we must subdivide it by the H=0 curves and parametrize separately each component of the complement. To follow a motion of the surface, one has to track the boundary shifts, mergers and splits of these components. On the other hand, Ushakov proposes in [113] an alternative, PDEbased, coordinate scheme viewing developable surfaces as solutions of the trivial Monge-Ampère equation. From an analytic viewpoint, developable surfaces can be seen locally, once they have been parametrized in the form z=z(x,y) , as solutions of the so-called trivial MongeAmpère partial differential equation, which is indeed the simplest of Monge-Ampère equations  zxx zxy zxy zyy  =0. (6.1) This is the case since a surface described as the graph of a function z=z(x,y), has Gaussian curvature: K=zxx ·zyy −z2 xy 1+z2 x+z2 y2. The PDE (6.1) almost provides a first set of coordinates for the space of torsal states U . If we have a developable surface S which is stably torsal, isometric to R , and admits a (non-isometric) parametrization by orthogonal projection to a plane then S admits a parametrization of the form z=z(x,y)and we can use 122 developable surfaces Theorem 3(Ushakov, [113]).The general solution to the equation (6.1) containing no flat patches is given in parametric form by x(u,v) = g(u)−v·f′(u)(6.2) y(u,v) = v(6.3) z(u,v) = u·g(u)−Zu 0g(t)dt +v· { f(u)−u·f′(u)}(6.4) where f(u) and g(u) are arbitrary functions such that f∈ C2,g∈ C1 , and g′(u)=0 everywhere. The functions f,g in Ushakov’s theorem give us a curve γ(u) = (g(u), 0, u·g(u)−Ru 0g(t)dt) and a vector field w(u) = (−f′(u), 1, f(u)− u·f′(u)) such that the ruled surface ϕ(u,v) = γ(u) + v·w(u) is actually developable. We could use (f,g) as coordinates for our stably torsal states, and explore the tangent space of these states as first-order deformations of the solutions but, unfortunately, the existence of a parametrization of the form z=z(x,y) can be assured only locally on a developable surface. Following the motion of the developable surface typically requires the dynamic subdivision of the original domain R in order to have such a parametrization, which leads even more intensely to the problem of tracking boundaries, mergers and splits of subdomains. 6.2 the boundary of developable surfaces There is an alternative approach to study theoretically the dynamics of developable surfaces isometric to a fixed bounded planar domain R : follow the motion of the boundary ∂R in space, and derive from this the developable surface that fills it. This leads to: Question. Given a piecewise smooth simple closed curve γ in R3 , what are the developable surfaces with boundary γ? The degeneracy nature of the trivial Monge-Ampère equation makes it fail to have a unique solution for this kind of boundary problem. Indeed, it is easy to find examples where there is more than one solution, as shown in Figure 6.1. Nevertheless, for problems such as the study of cloth dynamics, it is not necessary for the boundary problem to have a unique solution. It suffices to know that it will always have a finite set of solutions, because this solution set is then discrete, with different solutions separated by a nontrivial jump in any tagging energy, local coordinates, etc. In such case, once one has a developable ruling with a boundary γ0 at time t=0 , the evolution γt of the boundary will determine the analytic continuation of the t=0 developable ruling, and identify a unique ruling for every time t . Herein lies the interest of the authors in 6.2 the boundary of developable surfaces 123 Figure 6.1: A smooth simple closed curve (in black) which is the boundary of two developable surfaces (indicated in red and blue respectively) Theorem 4(Theorem).Let γ be a simple closed curve in R3 which is piecewise C2 , has nonvanishing curvature, its torsion vanishes at finitely many points, and such that only for finitely many pairs s=˜ s does the tangent line to γ at s pass through γ(˜ s) . Then, there can be at most finitely many developable surfaces with boundary γ and nonzero mean curvature in its interior. Let us point out that the preconditions that we impose on γ are generic, i.e. satisfied by a dense open subset of the embeddings of S1 in R3 . The starting idea to prove the theorem is another classical result, analogous to Proposition 2: Lemma 2.Let S be a torsal surface with boundary γ , and l⊂S a segment with endpoints P,Q in γ . Then the common tangent plane to S along l is tangent to γ at both P,Q , i.e. γ′(P) and γ′(Q) are both contained in the plane. The proof consists in pointing out that the normal vector to S stays constant over the segment l , and that γ is tangent to S . Lemma 2 presents developable rulings as arcs of bitangent planes (i.e., tangent to γ at 2points). Such planes are given by pairs s=˜ s whose tangent lines are coplanar: Proposition 4.Let γ:[0, L]⊂R3 be a simple, closed, arc-parametrized C3curve. The function D:[0, L]2−→ R (s,˜ s)7−→ det γ(s)−γ(˜ s),γ′(s),γ′(˜ s) is a Morse function at a neighborhood of its zeros (s,˜ s) such that: s=˜ s , γ has nonzero curvature and torsion at both s,˜ s , and the tangent line to γin each one does not pass through the other point of the curve. Proof. It is a straightforward computation. With coordinates (s,˜ s) we have that dD =det γ(s)−γ(˜ s),γ′′(s),γ′(˜ s), det γ(s)−γ(˜ s),γ′(s),γ′′(˜ s) 124 developable surfaces Let (s,˜ s) be a zero of D with s=˜ s , which is also a critical point of D . If any of the linear subspaces spanned by γ(s)−γ(˜ s),γ′(˜ s) and by γ(s)−γ(˜ s),γ′(s) has dimension less than 2, the tangent line to γ at one of the points γ(s),γ(˜ s) contains the other. When both linear subspaces have dimension 2, the conditions D(s,˜ s) = 0, dD(s,˜ s)=(0, 0) show that γ has the same osculating plane to γ at the points γ(s),γ(˜ s) . Because of this, the second differential d2Dof Dis κsτsdet (γ(s)−γ(˜ s),Bs,γ′(˜ s))0 0κ˜ sτ˜ sdet (γ(s)−γ(˜ s),γ′(s),B˜ s)! Here κ,τ,B are respectively the curvature, torsion, binormal vector of the Frenet frame, at the point given by their subindex. The determinants in the diagonal of d2D are nonzero because each consists of a binormal vector and a basis for the osculating plane at the same point of the curve. Proposition 4has a version for piecewise C3 curves, saying just that D is Morse under the additional hypothesis that s,˜ s do not correspond to corner points, at which D has two different definitions. We are now ready for Proof of Theorem 4. A torsal developable surface S is foliated by segments that can only end at the boundary or at points of vanishing mean curvature. Having ruled out the latter, S is determined by an arc of bitangent planes B(t) , with t∈[a,b] , such that the curves s(B(t)),˜ s(B(t)) formed by the two points of tangency of B(t) cover γ . Away from the finite set of horizontal and vertical lines in [0, L]2 where one of the values s,˜ s corresponds to a corner point, point with vanishing torsion, or point whose tangent line intersects γ again, the pairs of values s=˜ s for which there exists at all a bitangent plane to γ through γ(s),γ(˜ s) lie by Proposition 4in the zero set of a Morse function D from an open subset of [0, L]2⊂R2 to R . The function D is proper, therefore it is Morse over a suitably small range of values (−ε,ε) , which implies that D0 is a finite union of smooth curves with transverse intersections in [0, L]2 . The arc B(t) is determined by its tangency points curve (s(B(t)),˜ s(B(t))) ⊂[0, L]2, which must lie in the union of D0 and finitely many vertical and horizontal lines, and cover γ , i.e. γ=s(B)∪˜ s(B) . There are only finitely many possibilities for that, once we specify a beginning point for the curves s(B),˜ s(B). 7 SEMANTIC CLASSIFICATION OF CLOTH STATES One of the reasons why robotic manipulation of cloth is a very challenging task is the infinite-dimensional shape-state space of cloth, which makes its state estimation very difficult. In this chapter, we introduce the dGLI Cloth Coordinates, a low finite-dimensional representation of the state of a piece of cloth that allows us to efficiently distinguish key topological changes in a folding sequence. Our representation is based on a directional derivative of the well-known Gauss Linking Integral. The proposed dGLI Cloth Coordinates are shown to be more accurate in the separation of cloth states than other classic shape distance methods when applied to a rectangular cloth. In order to test this representation (and others), we use the full working simulator developed in the previous chapters to generate a full data-base of folded cloth states. After reviewing the state of the art in Section 7.1, we present preliminary concepts, such as the Gauss Linking Integral, and explain its limitations in a planar setting in Section 7.2. Then, in Section 7.3 we introduce the novel concept of the directional derivative of the GLI which is also applicable in a flat configuration. We derive first an expression for the GLI of two segments, then we prove that we can perturb the segments slightly to obtain information when they are co-planar and we explain how to apply this to a fully meshed cloth. Finally, in Section 7.4we apply this new index to a data-base of cloth configurations of a napkin taken from simulated folding sequences and then we test experimentally the index on real images of folded clothes. 7.1 related work As stated in the introduction, textile objects are important and omnipresent in many relevant scenarios of our daily lives, like domestic, healthcare, or industrial contexts. However, as opposed to rigid objects, whose pose is fixed with position and orientation, textile objects are challenging to handle for robots because they change shape under contact and motion, resulting in an infinite-dimensional configuration space. This huge dimensional jump makes existing perception and manipulation methods difficult to apply to textiles. Recent reviews on cloth manipulation, like [101,122], agree on the need to find a simplified representation that enables more powerful learning methods to solve different problems related to cloth manipulation. 125 126 semantic classification of cloth states Different representations have been used in the literature of cloth manipulation, e.g. silhouette representations [80] or contours [31], assuming the high-level reasoning on cloth states was given. More modern end-to-end learning approaches use RGB-D images as direct input [58,72,78,103,109], but only very simple actions can be defined due to the limited state representation. In addition, these methods need large amounts of real or simulated data that are expensive to obtain and label, as no underlying previous knowledge is used to understand the geometric relationship between different states. Therefore, finding a low-dimensional representation of cloth based on low-level features remains an active open problem, while the high-level aspect of understanding cloth deformation is still almost unexplored. Furthermore, to enable reasoning, abstraction and planning, rigid object manipulation applies object recognition methods in order to link objects to actions/affordances [16,119]. Contacts are estimated among the objects to recognize states such as “on top of", “inside of" [4]. However, when it comes to cloth manipulation, no work has explored the semantic state identification that could lead to particular actions depending on the task in mind. For simpler deformable objects like a box with an articulated lid, the open configuration clearly allows the action of closing the box or picking something from inside. An equivalent example for cloth would be to recognize a folded corner that needs to be either flattened back if the task is to lay it flat on the table, or picked up if the task is folding. In this context, we wish to classify the configuration space of a piece of cloth in macro-states (or just states), where each state is the set of cloth configurations that can be manipulated in the same way, i.e., that have similar grasping affordances. Figure 7.1: Folding sequence of a quadrangular cloth with its associated dGLI cloth coordinates, represented as upper triangular matrices. Each matrix element mij is a geometrical value corresponding to the dGLI between the segments i and j highlighted in red in the corresponding folded state. Notice how some values of the matrix change sign when corners are folded or cross each other. In this work, we present a coordinate representation of the configuration of a rectangular cloth as an upper triangular matrix form (see Figure 7.1). This representation can be computed with a closed-form 7.2 preliminaries 127 formula from low-level features of the cloth, mainly the position of its border, and enables the recognition and classification of high-level states, since we can define a distance between cloth configurations. That allows us to classify different configurations into states that we identify as different, meaning that they afford different actions. Our coordinates are based on a topological index, the Gauss Linking Integral ( GLI ). This index has been used in the past for robotic manipulation [57,91,107,108,124] but can only be applied to 3D curves. For a pair of almost coplanar curves, as the boundary curves of a folded garment, the GLI vanishes and it ceases to be informative. In order to consistently consider 2D curves as well as 3D curves, we introduce in this work the concept of the directional derivative of the GLI , dGLI , applied to a pair of curves. The dGLI is symmetric on the curves and it only depends on the relative position between them. We assign the dGLI Cloth Coordinates to a state of a garment as follows: first select a subset of edges (it may contain all of them) from a discretization of the boundary of the garment; then fix an ordering on these edges and compute the dGLI between any pair of edges in their spatial position of the current state of the garment; this gives a symmetric matrix from which only the upper triangular part is taken in order to avoid redundancies; the dGLI cloth coordinates of the state are precisely the entries of this upper triangular matrix (see Figure 7.1). Our resulting representation can be computed efficiently and is invariant under isometric movements of the garment (i.e. rotations and translations), leaving invariant a distinguished direction normal to a predominant plane in the scene (e.g. a table used as support for the manipulations). 7.2 preliminaries Given two non-intersecting 3D-space curves γ1 , γ2 parameterized by x(s) and y(t) , respectively, with s,t∈I= [0, 1] , the Gauss Linking Integral between them, GLI for short, is GLI(γ1,γ2) = 1 4πZIZI (y(t)−x(s)) ·(y′(t)×x′(s)) |y(t)−x(s)|3dtds or written in a compact way GLI(γ1,γ2) = 1 4πZ Z (γ2−γ1)·[γ′ 2×γ′ 1] ∥γ2−γ1∥3. (7.1) This double integral is invariant under re-parameterizations of the curves. In the case that both curves γ1 and γ2 are closed and smooth, their GLI is integer valued (due to the chosen normalization factor 1 4π ) and it is an invariant of the topology of the embedded curves (see [2]). Historically, the GLI was first introduced by Gauss, presumably related to his works on magnetism (according to [98]) or on astronomy 128 semantic classification of cloth states (according to [33]). Considering the GLI(γ,γ) of twice the same nonself-intersecting smooth curve γ , then the double integral (taking the domain of integration outside the diagonal of I×I ) defines another geometric invariant of the curve, known as writhe or writhing number of γ . Despite their resemblance, the GLI and the writhe measure different quantities: consider a normal vector field v of length ϵ>0 on γ , and the curve γv of endpoints of the vector field v , which is embedded and in one-to-one smooth correspondence with γ for sufficiently small ϵ . Then the GLI of these two close copies of the same γ differs from the writhe in GLI(γ,γv)−GLI(γ,γ) equal to the total twist of v . This result is known as the C˘alug˘areanu-White-Fuller theorem (see [90]). However, both indexes, GLI and writhe, are non-informative for planar curves, since they both vanish. The GLI has been used for many applications after a version of the above formula for polygonal curves appeared in the context of DNA protein structures [68], with additional efficient formulations given in [66], from which we have chosen the following: given a discretization of the curves into N and M segments, that is, γ1= {γPiPi+1,i=1, . . . , N} and γ2={γQiQi+1,i=1, . . . , M} , where each segment is parameterized as γAB(s) = A+s AB for s∈[0, 1] , then the GLI between both curves is GLI(γ1,γ2) = 1 4π N ∑ i=1 M ∑ j=1 GLI(γPiPi+1,γQiQi+1)(7.2) where the GLI between a pair of segments γAB and γCD is computed as GLI(γAB,γCD) = arcsin( nA nD) + arcsin( nD nB) +arcsin( nB nC) + arcsin( nC nA)(7.3) with  nA=∥ AC × AD∥, nB=∥ BD × BC∥,  nC=∥ BC × AC∥, and  nD=∥ AD × BD∥. The above discrete formula was used by Ho [54] to identify and synthesize animated characters in intertwined positions [53,54]. In the context of robotics, the GLI has been applied to representative curves of the workspace to guide path planning through holes [57,124], for guiding caging grasps in [91,107,108], and for planning humanoid robot motions, using the GLI to guide reinforcement learning [123]. In this work, for the first time, we develop a further analysis of the notion to be able to apply it to planar or almost planar curves, which opens the door to a wider spectrum of applications. 7.3 derivation of the cloth coordinates 129 7.3 derivation of the cloth coordinates As we have mentioned above, the GLI of two coplanar curves vanishes; so for many configurations of robotic interest — configurations where the cloth is nearly flat on a table, ready to be folded or already folded— the GLI does not provide much information. Our aim in this section is therefore to develop a similar index able to distinguish planar configurations. We shall see that a natural index to consider is in fact a directional derivative of the GLI , but to arrive at such an index we must first make a few observations about the GLI when applied to pairs of segments as in Equation (7.3); since the class of curves we will be working with computationally are piece-wise linear. 7.3.1GLI of two segments Since two segments AB and CD are uniquely defined by the four endpoints A,B,C,D∈R3 , the GLI of two segments computed in Equation (7.3) can be viewed as a function from (R3)4≡R12 to R . To emphasize that from now on we are considering segments we define G:R12 →Ras G(A,B,C,D) = GLI(γAB,γCD) = 1 4πZ Z (γCD −γAB)·[γ′ CD ×γ′ AB] ∥γCD −γAB∥3. (7.4) Note that technically G is not defined in the whole of R12 , since it is not defined when γAB and γCD intersect. Next, we will find a reformulation for Gwherever it is defined. Notice that the numerator in the integral expression of the GLI of Equation (7.4) is constant (for any tand s) and equals (γCD −γAB)·[γ′ CD ×γ′ AB] = = (  AC +t CD −s AB)·[ CD × AB] = = AC ·[ CD × AB] = = AC ·[(  CA + AD)× AB] = = AC ·[ AD × AB] =  AB ·[ AC × AD] = =det( AB, AC, AD) the signed volume of the tetrahedron ABCD multiplied by 6 . By writing V(A,B,C,D) = det( AB, AC, AD) and I(A,B,C,D) = 1 4πZ Z 1 ∥γCD −γAB∥3, we have G=V · I . (7.5)