Realtime Simulation of Thin-Shell Deformable Materials using CNN-Based Mesh Embedding
Qingyang Tan Zherong Pan Lin Gao Dinesh Manocha



  




















[1]▹ #1
SE[DOWHILE]DodoWhile[1] #1


#1#
  #1





FigFigure <ref>
figFigure <ref>
parSection <ref>
appenAppendix <ref>
secSection <ref>
subSection <ref>
tableTable <ref>
algAlgorithm <ref>
AlgAlgorithm <ref>
DefDefinition <ref>
ThmTheorem <ref>
stepStep <ref>
lnLine <ref>
eqEquation <ref>
pbProblem <ref>
itItem <ref>
teTerm <ref>
Eq:#1:(<ref>)
EqEquation #1:











































eqntiny
    eqntinyN
        @font5.74pt5.74pt1ruled@notopfs@cfontfs@capt@ruled
          fs@prefs@post2ptfs@mid2pt2ptfs@iftopcaptalgorithm@ruled@notop
         Realtime Simulation of Thin-Shell Deformable Materials 
        using CNN-Based Mesh Embedding
        
        
            Qingyang Tan^1, Zherong Pan^2, Lin Gao^3, and Dinesh Manocha^1
        
        Video Link: <https://youtu.be/zuXoQYJeAfc>
        ^1Qingyang Tan and Dinesh Manocha are with Department of Computer Science and Electrical & Computer Engineering, University of Maryland at College Park. {qytan,dm@cs.umd.edu}^2Zherong Pan is with Department of Computer Science, University of North Carolina at Chapel Hill. {zherong@cs.unc.edu}^3Lin Gao is with Institute of Computing Technology, Chinese Academy of Sciences. {gaolin@ict.ac.cn}
            December 30, 2023
        ========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================
        
        emptyempty
        We address the problem of accelerating thin-shell deformable object simulations by dimension reduction. We present a new algorithm to embed a high-dimensional configuration space of deformable objects in a low-dimensional feature space, where the configurations of objects and feature points have approximate one-to-one mapping. Our key technique is a graph-based convolutional neural network (CNN) defined on meshes with arbitrary topologies and a new mesh embedding approach based on physics-inspired loss term. We have applied our approach to accelerate high-resolution thin shell simulations corresponding to cloth-like materials, where the configuration space has tens of thousands of degrees of freedom. We show that our physics-inspired embedding approach leads to higher accuracy compared with prior mesh embedding methods. Finally, we show that the temporal evolution of the mesh in the feature space can also be learned using a recurrent neural network (RNN) leading to fully learnable physics simulators. After training our learned simulator runs 500-10000× faster and the accuracy is high enough for robot manipulation tasks.
        
        
        § INTRODUCTION
        
        A key component in robot manipulation tasks is a dynamic model of target objects to be manipulated. Typical applications include cloth manipulation <cit.>, liquid manipulation <cit.>, and in-hand rigid object manipulation <cit.>. Of these objects, cloth is unique in that it is modeled as a thin-shell, i.e., a 2D deformable object embedded in a 3D workspace. To model the dynamic behaviors of thin-shell deformable objects, people typically use high-resolution meshes (e.g. with thousands of vertices) to represent the deformable objects. Many techniques have been developed to derive a dynamic model under a mesh-based representation, including the finite-element method <cit.>, the mass-spring system <cit.>, the thin-shell model <cit.>, etc. However, the complexity of these techniques can vary from 𝒪(n^1.5) to 𝒪(n^3)<cit.>, where n is the number of DOFs, which makes them very computationally cost on high-resolution meshes. For example, <cit.> reported an average computational time of over 1 minute for predicting a single future state of a thin-shell mesh with around 5000 vertices. This simulation overhead is a major cost in various cloth manipulation algorithms including <cit.>.
        
        In order to reduce the computational cost, one recent trend is to develop machine learning methods to compute low-dimensional embeddings of these meshes. Low-dimensional embeddings were original developed for applications such as image compression <cit.> and dimension reduction <cit.>. The key idea is to find a low-dimensional feature space with approximate one-to-one mapping between a low-dimensional feature point and a high-dimensional mesh shape. So that the low-dimensional feature point can be treated as an efficient, surrogate representation of the original mesh. 
        
        However, computing low-dimensional embeddings for general meshes poses new challenges because, unlike 2D images, meshes are represented by a set of unstructured vertices connected by edges and these vertices can undergo large distortions when cloth deforms. As a result, a central problem in representing mesh deformation data is to find an effective parameterization of the feature space that can handle arbitrary mesh topologies and large, nonlinear deformations. Several methods for low-dimensional mesh embeddings are based on PCA <cit.>, localized PCA <cit.>, and Gaussian Process <cit.>. However, these methods are based on vertex-position features and cannot handle large deformations.
        
        Main Results: We present a novel approach that uses physics-based constraints to improve the accuracy of low-dimensional embedding of arbitrary meshes for deformable simulation. We further present a fully learnable physics simulator of clothes in the feature space. The novel components of our algorithm include:
        
          *  A mesh embedding approach that takes into account the inertial and internal potential forces used by a physical simulator, which is achieved by introducing a physics-inspired loss function term, i.e., vertex-level physics-based loss term (PB-loss). This also preserves the material properties of the mesh.
        
          *  A stateful, recurrent feature-space physics simulator that predicts the temporal changes of meshes in the feature space, which are modeled by introducing accurate enough for learning cloth features and training cloth manipulation controllers (see fig:simulateRobot).
        
        To test the accuracy of our method, we construct multiple datasets by running cloth simulations using a high-resolution mesh under different material models, material parameters, and mesh topologies. We show that our embedding approach leads to better accuracy in terms of physics rule preservation than prior method <cit.> that uses only a data term, with up to 70% improvement. We have also observed up to 19% and 18% improvements in mesh embedding accuracy on commonly used metrics such as ℳ_rms and ℳ_STED. Finally, we show that our feature space physics simulator can robustly predict dynamic behaviors of clothes undergoing unseen robot manipulations, while achieving 500-10000× speedup over simulators running in the high-dimensional configuration space.
        
        
        The paper is organized as follows. We first review related work in sec:related. We define our problem and introduce the basic method of low-dimensional mesh embedding in sec:VDE. We introduce our novel PB-loss and the learnable simulation architecture in sec:physics. Finally, we describe the applications in sec:app and highlight the results in sec:results.
        
        
        § RELATED WORK AND BACKGROUND
        
        We summarize related work in mesh deformations and representations, deformable object simulations, and machine learning methods for mesh deformations. 
        
        Deformable Simulation for Robotics are frequently encountered in service robots applications such as laundry cleaning <cit.> and automatic cloth dressing <cit.>. Studying these objects can also benefit the design of soft robots <cit.>. While these soft robots are usually 3D volumetric deformable objects, we focus on 2D shell-like deformable objects or clothes. In some applications such as visual servoing <cit.> and tracking <cit.>, deformable objects are represented using point clouds. In other applications including model-based control <cit.> and reconstruction <cit.>, the deformable objects are represented using meshes and their dynamics are modeled by discretizing the governing equations using the finite element method (FEM). Solving the discretized governing equation is a major bottleneck in training a cloth manipulation controller, e.g., <cit.> reported up to 5 hours of CPU time spend on thin-shell simulation which is 4-5 times more costly than the control algorithm.
        
        Deformable Object Simulations is a key component in various model-based control algorithms such as virtual surgery <cit.> and soft robot controllers <cit.>. However, physics simulators based on the finite element method <cit.>, the boundary-element method <cit.>, or simplified models such as the mass-spring system <cit.> have a superlinear complexity. An analysis is given in <cit.>, resulting in 𝒪(n^1.5) complexity, where n is the number of DOFs. In a high-resolution simulation, n can be in the tens of thousands. As a result, learning-based methods have recently been used to accelerate physics simulations. This can be done by simulating under a low-resolution using FEM and then upsampling <cit.> or by learning the dynamics behaviors of clothes <cit.> and fluids <cit.>. However, these methods are either not based on meshes <cit.> or not able to handle arbitrary topologies <cit.>.
        
        Machine Learning Methods for Mesh Deformations has been in use for over two decades, of which most methods are essentially low-dimensional embedding techniques. Early work are based on principle component analysis (PCA) <cit.> that can only represent small, local deformations or Gaussian processes <cit.> that are computationally costly to train and do not scale to large datasets. Recently, deep neural networks have been used to embed high-dimensional nonlinear functions <cit.>. However, these methods rely on regular data structures such as 2D images. To handle meshes with arbitrary topologies, earlier methods <cit.> represent a mesh as a 3D voxelized grid or reconstruct 3D shapes from 2D images <cit.> using a projection layer. Recently, methods have been proposed to define CNN directly on mesh surfaces, such as CNN on parametrized texture space <cit.>, and CNN based on spatial filtering <cit.>. The later has been used in <cit.> to embed large-scale deformations of general meshes. Our contribution is orthogonal to these techniques and can be used to improve the embedding accuracy for any one of these methods.
        
        
        § LOW-DIMENSIONAL MESH EMBEDDING
        
        In this section, we provide an overview of low-dimensional embedding of thin shell like meshes such as clothes. Our goal is to represent a set of N deformed meshes, S_m, with each mesh represented using a set of K vertices, denoted as p_m∈R^3K. We denote the ith vertex as p_m,i∈R^3. Here m=1,⋯,N and i=1,⋯,K. These vertices are connected by edges, so we can define the 1-ring neighbor set, 𝒩^1_i, and the 2-ring neighbor set, 𝒩^2_i, for each p_i, as shown in fig:method (c). Our goal is to find a map z→p, where z is a low-dimensional feature and p∈R^3K such that, for each m, there exists a z_m where z_m is mapped to a mesh close to p_m. To define such a function, we use graph-based CNN and ACAP features <cit.> to represent large-scale deformations.
        
        
        
         §.§ ACAP Feature
        
        For each S_m, an ACAP feature is computed by first finding the deformation gradient T_m,i on each vertex:
        
            .4!T_m,i≜Tj∈𝒩^1_ic_ij(p_m,i-p_m,j) - T(p_1,i-p_1,j)^2,
        
        where c_ij are cotangent weights <cit.>. Here, we use S_1 as a reference shape. Next, we perform polar decomposition to compute T_m,i=R_m,iS_m,i where R_m,i is orthogonal and S_m,i is symmetric. Finally, R_m,i is transformed into log-space in an as-consistent-as-possible manner using mixed-integer programming. The final ACAP feature is defined as ACAP_m,i≜{log(R_m,i),S_m,i}∈R^9 due to the symmetry of S_m,i. We denote the ACAP feature transform as: ACAP(p_m)∈R^9K. It is suggested, e.g., in <cit.>, that mapping z_m to the ACAP feature space leads to better effectiveness in representing large-scale deformations. Therefore, we define our mapping function to be D(z): z→ACAP(p) and then recover p via the inverse feature transform: ACAP^-1.
        
        
        
         §.§ Graph-Based CNN for Feature Embedding
        
        The key idea in handling arbitrary mesh topologies is to define D as a graph-based CNN using local filters <cit.>:
        
            D≜C_L^T⋯C_2^T∘C_1^T∘F^T,
        
        where L is the number of convolutional layers and C^T is the transpose of a graph-based convolutional operator. Finally, F is a fully connected layer. Each layer is appended by a leaky ReLU activation layer. A graph-based convolutional layer is a linear operator defined as:
        
            .45!C(ACAP(p_m))_m,i≜W×ACAP_m,i+
            W_𝒩j∈𝒩^1_iACAP_m,j/|𝒩^1_i|+b,
        
        where W,W_𝒩,b are optimizable weights and biases, respectively. All the weights in the CNN are trained in a self-supervised manner using an autoencoder and the reconstruction loss:
        
            .42!ℒ_recon=m=1ND∘E∘ACAP(p_m)-ACAP(p_m)^2/N,
        where E is a mirrored encoder of D with a weight-tied architecture defined as:
        
            E≜F∘C_1⋯C_L-1∘C_L,
        
        which means that each layer in E is a transpose of the corresponding layer in D with shared weights. The construction of this CNN is illustrated in fig:method (a). In the next section, we extend this framework to make it aware of physics rules.
        
        
        § PHYSICS-BASED LOSS TERM
        
        We present a novel physics-inspired loss term that improves the accuracy of low-dimensional mesh embedding. Our goal is to combine physics-based constraints with graph-based CNNs, where our physics-based constraints take a general form and can be used with any material models such as FEM <cit.> and mass-spring system <cit.>. We assumes that S_m is generated using a physics simulator that solves a continuous-time PDE of the form:
        
            Mpt=-force(p,q),
        
        where M is the mass matrix and t is the time. This form of governing equation is the basis for state-of-the-art thin shell simulators including <cit.>. force(p,q) models internal and external forces affecting the current mesh p. The force is also a function of the current control parameters q, which are the positions of the grasping points on the mesh (red dots of fig:datasetVis). This continuous time PDE eq:PDE can be discretized into N timesteps such that S_m is the position of S at time instance iΔ t, where Δ t is the timestep size. A discrete physics simulator can determine all p_m given the initial condition p_1,p_2 and the sequence of control parameters q_1,q_2,⋯,q_N by the recurrent function:
        
            p_m≜ f(p_m-2,p_m-1,q_m),
        
        where f is a discretization of eq:PDE. To define this discretization, we use a derivation of <cit.> that reformulates f as the following optimization:
        
            p_m   ≜pℒ_phys(p_m-2,p_m-1,p,q_m)
        
            ℒ_phys   ≜p-2p_m-1+p_m-2_M^2/(2Δ t^2)+𝒫(p,q_m).
        
        Note that eq:PLOSS is just one possible implementation of eq:SIM. Here the first term models the kinematic energy, which requires each vertex to move in its own velocity as much as possible if no external forces are exerted. The second term models forces caused by various potential energies at configuration p. In this work, we consider three kinds of potential energy: 
        
          *  Gravitational energy 𝒫_g(p)≜-i=1Kg^TMp, where g is the gravitational acceleration vector.
            
          *  Stretch resistance energy, 𝒫_s, models the potential force induced by stretching the material.
            
          *  Bending resistance energy, 𝒫_b, models the potential force induced by bending the material.
        
        There are many ways to discretize 𝒫_s,𝒫_b, such as the finite element method used in <cit.> or the mass-spring model used in <cit.>. Both formulations are evaluated in this work. 
        
          * <cit.> models the stretch resistance term, 𝒫_s, as a set of Hooke's springs between each vertex and vertices in its 1-ring neighbors. In addition, the bend resistance term, 𝒫_b, is defined as another set of Hooke's springs between each vertex and vertices in its 2-ring neighbors. (fig:method (c))
            
          * <cit.> models the stretch resistance term, 𝒫_s, as a linear elastic energy resisting the in-plane deformations of each mesh triangle. In addition, the bend resistance term, 𝒫_b, is defined as a quadratic penalty term resisting the change of the dihedral angle between any pair of two neighboring triangles. (fig:method (d))
        
        
        Our approach uses eq:PLOSS as an additional loss function for training D,E. Since eq:PLOSS is used for data generation, using it for mesh deformation embedding should improve the accuracy of the embedded shapes. However, there are two inherent difficulties in using 𝒫 as an loss function. First, 𝒫 is defined on the vertex level as a function of p_m, not on the feature level as a function of ACAP(p_m). To address this issue, we use the inverse function ACAP^-1 to reconstruct p_m from ACAP(p_m). The implementation of ACAP^-1 is introduced in sec:IssueA. By combining ACAP^-1 with ℒ_phys, we can train the mesh deformation embedding network using the following loss:
        
            ℒ̃_phys≜ℒ_phys(ACAP^-1∘D(z_m, m-1, m-2), q_m)
               ℒ_ephys=ℒ̃_phys(E∘ACAP(p_m,m-1,m-2),q_m).
        
        Our second difficulty is that the embedding network is stateless and does not account for temporal information. In other words, function E only takes p_m as input, while eq:PLOSS requires p_m,p_m-1,p_m-2. To address this issue, we use a small, fully connected, recurrent network to represent the physics simulation procedure in the feature space. The training of this stateful network is introduced in sec:IssueB. Finally, in addition to the PB-loss, we also add an autoencoder reconstruction loss on the vertex level as a regularization:
        
            .45!ℒ_vert=m=1NACAP^-1∘D∘E∘ACAP(p_m)-p_m^2/N.
        This vertex level loss can be removed from our loss function without significantly affect the quality of results. However, ℒ_vert provides complementary information to ℒ_recon. Since ACAP is feature in the gradient domain, using only ℒ_recon will reconstruct accurate local geometric features, but can lead to large error in vertices' positions. Therefore, we combine ℒ_vert and ℒ_recon to reduce the gradient domain errors and absolute vertex position errors.
        
         §.§ The Inverse of the ACAP Feature Extractor
        
        The inverse of the ACAP function (black block in fig:method) involves three steps. Fortunately, each step can be easily implemented in a modern neural network toolbox such as TensorFlow <cit.>. The first step computes R_m,i from Log(R_m,i) using the Rodrigues' rotation formula, which involves only basic mathematical functions such as dot-product, cross-product, and the cosine function. The second step computes T_m,i from R_m,i,S_m,i, which is a matrix-matrix product. The final step computes p_m,i from T_m,i. According to eq:RECON, this amounts to pre-multiplying the inverse of a fixed sparse matrix, L, representing the Poisson reconstruction. However, this L is rank-3 deficient because it is invariant to rigid translation. Therefore, we choose to define a pseudo-inverse by fixing the position of the grasping points q:
        
            L^†p≜I0LA^TA0^-1pq,
        
        which can be pre-factorized. Here A^3×3K is a matrix selecting the grasping points.
        
        
        
         §.§ Stateful Recurrent Neural Network
        
        A physics simulation procedure is Markovian, i.e. current configuration p_m-1p_m only depends on previous configuration p_m-2p_m-1 of the mesh. As a result, ℒ_phys is a function of both p_m-2, p_m-1, and p_m, which measures the violation of physical rules. However, our embedding network is stateless and only models p_m. In order to learn the entire dynamic behavior, we augment the embedding network with a stateful, recurrent network represented as a multilayer perceptron (MLP). This MLP represents a physically correct simulation trajectory in the feature space and is also Markovian, denoted as:
        
            MLP(z_m-2,z_m-1,q_m)=z_m.
        
        Here the additional control parameters q are given to MLP as additional information. 
        We can build a simple reconstruction loss below to optimize MLP:
        
            .45!ℒ_sim = m=3NMLP(z_m-2,z_m-1,q_m)-z_m^2/(N-2).
        
        In addition, we can also add PB-loss to train this MLP, for which we define ℒ_mphys on a sequence of N meshes by unrolling the recurrent network:
        
            ℒ_mphys = 1/N-2* 
            
            (   ℒ̃_phys(z_1,z_2,MLP(z_1,z_2,q_3),q_3)+
               ℒ̃_phys(z_2,z_3,MLP(z_2,z_3,q_4),q_4)+⋯+
               ℒ̃_phys(z_N-2,z_N-1,MLP(z_N-2,z_N-1,q_N),q_N)).
        
        However, we argue that eq:LPHYS will lead to a physically incorrect result and cannot be directly used for training. To see this, we note that eq:PLOSS is the variational form of eq:PDE. So that p_m is physically correct when ℒ_phys is at its local minima, i.e. the following partial derivative vanishes:
        
            ℒ_phys(p_m-2,p_m-1,p_m,q_m)p_m=0  ∀ m.
         
        However, if we sum up ℒ_phys over a sequence of N meshes and require the summed-up loss to be at a local minimum, as is done in eq:LPHYS, then we are essentially requiring the following derivatives to vanish:
        
            ℒ_phys(p_m-2,p_m-1,p_m,q_m)p_m+
               ℒ_phys(p_m-1,p_m,p_m+1,q_m)p_m+
               ℒ_phys(p_m,p_m+1,p_m+2,q_m)p_m=0  ∀ m.
        
        
        The difference between eq:VIO_MEASURE and eq:VIO_MEASURE_SUM is the reason that eq:LPHYS gives an incorrect result. To resolve the problem, we slightly modify the back propagation procedure of our training process by setting the partial derivatives of ℒ_phys with respect to its first two parameters to zero:
        
            ℒ_phys(p_m-2,p_m-1,p_m,q_m)[p_m-1,p_m-2]=0  ∀ m,
        
        which, combined with eq:VIO_MEASURE_SUM, leads to eq:VIO_MEASURE. (We add similar gradient constraints when optimizing over eq:embedding_phys.) This procedure is equivalent to an alternating optimization procedure, where we first compute a sequence of feature space coordinates, z_m, using the recurrent network (eq:MLP) and then fix the first two parameters z_m-2,z_m-1 and optimize ℒ_mphys with respect to its third parameter z_m.
        
        
        
        
        § APPLICATIONS
        
        The two novel components in our method, the ACAP^-1 operator and the stateful PB-loss, enable a row of new applications, including realtime cloth inverse kinematics and feature space physics simulations. 
        
        
        
         §.§ Cloth Inverse Kinematics
        
        Our first application allows a robot to grasp several points of a piece of cloth and then infer the full kinematic configuration of the cloth. Such inverse kinematics can be achieved by minimizing a high-dimensional nonlinear potential energy, such as ARAP energy <cit.>, which is computationally costly. Using the inverse of the ACAP feature extractor, our method allows vertex-level constraints. Therefore, we can perform solve for the cloth configuration by a fast, low-dimensional minimization in the feature space as:
        
            z^*=zℒ∘ACAP^-1∘D(z),
         
        where we treat all the grasped vertices as control parameters q used in eq:invACAP. This application is stateless and the user controls a single feature of a mesh, z, so that we drop the kinetic term in ℒ_phys and only retain the potential term 𝒫_s+𝒫_b. Some inverse kinematic examples generated using this formulation are shown in fig:control. Note that detailed wrinkles and cloth-like deformations are synthesized in unconstrained parts of the meshes.
        
        
         §.§ Feature Space Physics Simulation
        
        For our second application, we approximate an entire cloth simulation sequence (eq:SIM) in the 128-dimensional feature space. Starting from z_1,z_2, we can generate an entire sequence of N frames by using the recurrent relationship in eq:MLP and can recover the meshes via the function ACAP^-1∘D. Such a latent space physics model has been previously proposed in <cit.> for voxelized grids, while our model works on surface meshes. We show two synthesized simulation sequences in fig:simulate. 
        
        
        
         §.§ Accuracy of Learned Simulator for Robotic Cloth Manipulation
        
        We show three benchmarks (fig:simulateRobot) from robot cloth manipulation tasks defined in prior work <cit.>. In these benchmarks, the robot is collaborating with human to maintain a target shape of a piece of cloth. To design such a collaborating robot controller, we use imitation learning by teaching the robot to recognized cloth shapes under various, uncertain human movements. Our learnable simulator can be used to efficiently generate these cloth shapes for training the controller. To this end, we train our neural-network using the original dataset from <cit.> obtained by running the FEM-based simulator <cit.>, which takes 3 hours. During test time, we perturb the human hands' grasp points along randomly directions. Our learned physical model can faithfully predict the dynamic movements of the cloth.
        
        
        § RESULTS
        
        To evaluate our method, we create two datasets of cloth simulations using eq:PLOSS. Our first dataset is called SHEET, which contains animations of a square-shaped cloth sheet swinging down under different conditions, as shown in fig:datasetVis (a). This dataset involves 6 simulation sequences, each with N=2400 frames. Among these 6 sequences, the first sequence uses the mass-spring model <cit.> to discretize eq:PDE and the cloth mesh has no holes (denoted as SHEET+<cit.>). The second and the third simulation sequences in the dataset use different material parameter, by multiplying the stretch/bend resistance term by 0.1 and thereby making the material softer and less resilient when stretched or bent. These two sequences are denoted as SHEET+<cit.>+0.1𝒫_s and SHEET+<cit.>+0.1𝒫_b, respectively.
        The forth sequence uses the mass spring model and the cloth mesh has holes, as shown in fig:simulate (a,b), which is denoted as (SHEET+<cit.>+holes). The fifth sequence uses FEM <cit.> to discretize eq:PDE and the cloth mesh has no holes (denoted as SHEET+<cit.>). The sixth sequence uses FEM to discretize eq:PDE and the cloth interacts with an obstacle, as shown in fig:datasetVis (c) (denoted as SHEET+<cit.>+obstacle). In the SHEET dataset, the cloth mesh without holes has K=4225 vertices and the cloth mesh with holes has K=4165 vertices. Our second dataset is called BALL, which contains animations of a cloth ball being dragged up and down under different conditions, as shown in fig:datasetVis (d). This dataset also involves 4 simulation sequences, each with N=500 frames. Using the same notation as the SHEET dataset, the 4 sequences in the BALL dataset are (BALL+<cit.>, BALL+<cit.>, BALL+<cit.>+0.1𝒫_s, BALL+<cit.>+0.1𝒫_b). 
        
        In the BALL dataset, the cloth ball mesh has K=1538 vertices. During comparison, for each dataset, we select first 12 frames in every 17 frames to form the training set. The other frames are used as the test set.
        
         §.§ Implementation
        
        We implement our method using Tensorflow <cit.> and we implement the PB-loss as a special network layer. When there is an obstacle interacting with the cloth, we model the collision between the cloth and the obstacle using a special potential term proposed in <cit.>. For better conditioning and a more robust initial guess, our training procedure is broken into three stages. During the first stage, we use the loss:
        
            ℒ_1=∑λ_iℒ_i  i∈{recon,vert,ephys}
        
        to optimize E,D. During the second stage, we use the loss:
        
            ℒ_2=∑λ_iℒ_i  i∈{sim,mphys}
        
        to optimize MLP. Finally, we add a fine-tuning step and use the loss: 
        
            ℒ_3=∑λ'_iℒ_i  i∈{recon,vert,ephys,sim}
        
        to optimize both E,D and MLP. Notice that, in order to train the mesh embedding network and MLP at the same time, we feed:
        
            z' = 0.5*z_m + 0.5*MLP(z_m-2, z_m-1, q_m)
        
        to D for better stability during the third stage.
        
        
        
         §.§ Physics Correctness of Low-Dimensional Embedding
        
        We first compare the quality of mesh deformation embeddings using two different methods. The quality of embedding is measured using three metrics. The first metric is the root mean square error, ℳ_rms<cit.>, which measures the averaged vertex-level error over all N shapes and K vertices. Our second metric is the STED metric, ℳ_STED<cit.>. This metric linearly combines several aspects of errors crucial to visual quality, including relative edge length changes and temporal smoothness. However, since ℳ_STED is only meaningful for consecutive frames, we compute ℳ_STED for the 5 consecutive frames in every 17 frames, which is the test set. Finally, we introduce a third metric, physics correctness, which measures how well the physics rule is preserved. Inspired by eq:VIO_MEASURE, physics correctness is measured by the norm of partial derivatives of ℒ_phys: ℳ_phys=ℒ_physp_m^2. Note that the absolute value of ℳ_phys can vary case by case. For example, ℳ_phys using the FEM method can be orders of magnitude larger than that using the mass-spring system in our dataset. So that only the relative value of ℳ_phys indicates improvement in physics correctness.
        
        
        Our first experiment compares the accuracy of mesh embedding with or without PB-loss. The version without PB-loss is our baseline, which is equivalent to adding vertex level loss to <cit.>. In addition, we remove the sparsity regularization from <cit.> to make it consistent with our formulation. We denote this baseline as <cit.>+ℒ_vert. A complete summary of our experimental results is given in table:comparisonA. The benefit of three-stage training is given in table:MLP. From table:comparisonA, we can see that including PB-loss significantly and consistently improves ℳ_phys. This improvement is large, up to 70% on the SHEET+<cit.> dataset. In addition, by adding ℒ_ephys, our method also better recognizes the relationship between each model and embeds them, thus improves ℳ_rms in all the cases. However, our method sometimes sacrifice ℳ_STED as temporal smoothness is not modeled explicitly in our method. Finally, we have added two rows to table:comparisonA, comparing our method with and without ℒ_vert, which shows that ℒ_vert effectively reduces the error in terms of absolute vertex positions.
        
         §.§ Discriminability of Feature Space
        
        In our second experiment, we evaluate the discriminability of mesh embedding by classifying the meshes using their feature space coordinates. Note that our datasets (fig:datasetVis) are generated by moving the grasping points back and forth. We use these movement directions as the labels for classification. For the SHEET dataset, we have 6 labels: ± X/Y,± R, where ± R means rotating the grasping points around ± Z axes. For the BALL dataset, we have 2 labels: ± Z. Note that it is trivial to classify the meshes if we know the velocity of the grasping points. However, this information is missing in our feature space coordinates because ACAP features are invariant to global rigid translation, which makes the classification challenging. 
        
        We visualize the feature space using t-SNE <cit.> compressed to 2 dimensions in fig:classify.
        We report retrieval performance in the KNN neighborhoods across different K’s, using method suggested by <cit.>. The normalized discounted cumulative gain (DCG) on the test set for SHEET+<cit.> is 0.8045 and for BALL+<cit.> is 0.9128.
        
        
        
        
         §.§ Sensitivity to Training Parameters
        
        In our third experiment, we evaluate the sensitivity of our method with respect to the weights of loss terms, as summarized in table:comparisonB. Our method outperforms <cit.>+ℒ_vert under a range of different parameters. We have also compared our method with other baselines such as <cit.> and <cit.>. As shown in the last two columns of table:comparisonB, they generate even worse result, which indicates that <cit.>+ℒ_vert is the best baseline.
        
        
         §.§ Robustness to Mesh Resolutions
        
        In this experiment, we highlight the robustness of our method to different mesh resolutions by lowering the resolution of our dataset. For SHEET+<cit.>, we create a mid-resolution counterpart with K=1089 vertices and a low-resolution counterpart with K=289 vertices. On these two new datasets, we compare the accuracy of mesh embedding with or without PB-loss. The results are given in table:comparisonC. Including PB-loss consistently improves ℳ_phys and overall embedding quality, no matter the resolution used.
        
        
         §.§ PB-Loss with Alternative Neural Network Architectures
        
        Our PB-loss is orthogonal to the architecture of neural networks. Therefore, we have conducted an additional set of experiments to highlight the performance improvement using a full-connected underlying neural network like <cit.>. The results are shown in table:fc.
        
         §.§ Difficulty in Contact Handling
        
        One exception appears in the SHEET+<cit.>+obstacle (blue row in table:comparisonA), where our method deteriorates physics correctness. This is the only dataset where the mesh is interacting with an obstacle. The deterioration is due to the additional loss term penalizing the penetration between the mesh and the obstacle. This term is non-smooth and has very high value and gradient when the mesh is in penetration, making the training procedure unstable. This means that direct learning a feature mapping for meshes with contacts and collisions can become unstable. However, we can solve this problem using a two-stage method, where we first learn a feature mapping for meshes without contacts and collisions, and then handle contacts and collisions at runtime using conventional method <cit.>, as is done in <cit.>.
        
        
        
         §.§ Speedup over FEM-Based Cloth Simulators
        In table:RUNNING_TIME, we have used our  neural network as a physics simulator and compared its performance with a conventional FEM-based method. Our method achieves 500-10000× speedup over the FEM-based method <cit.> on average.
        
        § CONCLUSION & LIMITATIONS
        
        In this paper, we present a new method that bridges the gap between mesh embedding and and physical simulation for efficient dynamic models of clothes. We achieve low-dimensional mesh embedding using a stateless, graph-based CNN that can handle arbitrary mesh topologies. To make the method aware of physics rules, we augment the embedding network with a stateful feature space simulator represented as a MLP. The learnable simulator is trained to minimize a physics-inspired loss term (PB-loss). This loss term is formulated on the vertex level and the transformation from the ACAP feature level to the vertex level is achieved using the inverse of the ACAP feature extractor.
        
        Our method can be used for several applications, including fast inverse kinematics of clothes and realtime feature space physics simulation. We have evaluated the accuracy and robustness of our method on two datasets of physics simulations with different material properties, mesh topologies, and collision configurations. Compared with previous models for embedding, our method achieves consistently better accuracy in terms of physics correctness and the mesh change smoothness metric (<cit.>).
        
        A future research direction is to apply our method to other kinds of deformable objects, i.e., volumetric objects <cit.>. Each and every step of our method can be trivially extended to handle volumetric objects by replacing the triangle surface mesh with a tetrahedral volume mesh. A minor limitation of the current method is that the stateful MLP and the stateless mesh embedding cannot be trained in a fully end-to-end fashion. We would like to explore new optimization methods to train the two networks in an end-to-end fashion while achieving good convergence behavior. Finally, our approach is limited to a single setting of thin-shell simulation and needs to be re-trained whenever there are changes in the resolution of the mesh, the material parameters, or the obstacles in the environment.
        
        § ACKNOWLEDGEMENT
        
        This work was supported by ARO grants (W911NF1810313 and  W911NF1910315) and Intel. Lin Gao was partially supported by the National Natural Science Foundation of China (No. 61872440) and Beijing Municipal Natural Science Foundation (No. L182016).
        
        IEEEtranS
