Handling deformable objects is a fundamental skill in robotic surgery, where surgeons precisely and minimally invasively treat soft tissue disease with robots. As the scope of robotic surgery expands, and as the gap between trained doctors and under-served patient populations widen, AIdriven surgical capabilities like autonomous suturing and tissue retraction become increasingly desirable and potentially necessary [1], [2]. Ultimately, these robots will require an understanding of tissue physics and deformable object manipulation in general in order to provide assistance or take over a surgical task.
Physics-based simulations have been proven to be a promising technique for deformable object manipulation (DOM) [3], [4], [5], [6]. Extensive research has been done on modeling, motion planning, and data representations for perception perspectives [7]. The development of physical simulation environments can support a variety of deformable objects, including thin-shell fabric [8], linear elastic ropes [9], volumetric tissue [10], [11], and fluids [12].
Considering many existing deformable simulators, the high computational cost and "reality gap" are two factors that Equal contributions. Authors are affiliated with Department of Electrical and Computer Engineering, University of California San Diego, La Jolla, CA 92093 USA.{x5liang, f4liu, yuz049, yul189, shl102, yip}@ucsd.edu It reveals suitable constraints' stiffness parameters for a PBD simulation of a deformable tissue in an online manner. Arrows in purple shows the immediate control direction at the given time step. The tissue's blue regions indicate lower stiffness, whereas red regions have higher stiffness. In this figure, the tissue is initialized with large stiffness parameters. limit the conventional simulation approaches. These include Finite Element Methods (FEM) and mass-spring system [13]. Uncertainties, inaccurately calibrated parameters, and unmodeled physical effects can all lead to the gap between simulation and reality. As a result, many simulators need their parameters fine-tuned online in order to be deployable for real-world robotic tasks. With the intention of bridging the gap between reality and simulation, the phrase real-to-sim was first used in [10] to describe approaches that compensate for the simulated errors from live observations.
Our goal is to investigate the real-to-sim gaps by taking into account both known physical (geometric, mechanical) constraints describing tissue deformation and observable high-dimensional data, i.e., point cloud. For real-time applications such as surgical tissue manipulation, grasping, and retraction, it would also be necessary to identify the proper simulation parameters of soft bodies with a fast online simulation approach.
The real-to-sim problem has recently gained attention in the literature. Most of the existing works have been focusing on closing the gap through effective policy transfer in a reinforcement learning (RL) setup, as reviewed in [14], [15]. Most of these works only rely on simulating rigid objects in the scene and robots with kinematics [16], [17]. Recently, several papers used similar deep RL for deformable objects, such as cloth [18], tissue [19], and ropes [20]. However, these works don't create explicit physical models but learn the system parameters or controls in an end-to-end manner.
TABLE I. SUMMARY OF SELECTED PREVIOUS WORKS THAT ADDRESSED DEFORMABLE REAL-TO-SIM GAP. Methods Residual Parameters Update Object types [10] 6 : volumetric [21] : 6 thin-shell [25] 6 : volumetric Proposed 6 6 both
It leads to several typical problems with learning strategies, including generalizability and data-hunger.
Other methods rely on physics-based simulation to minimize the real-to-sim gap. Table I presents a summary. In our previous paper [10], we directly update the simulated positions of the volumetric particles using the spatial gradient of the signed distance field of point cloud observation. However, the simulation parameters are not updated, necessitating frame-by-frame registration. Similarly, [21] optimizes for a simulation parameter, such as mass or stiffness, to minimize the difference between a simulation and the point cloud observation. Furthermore, the authors enhanced their approach in [22] as probabilistic inference over simulation parameters of the deformable object. However, their method is limited to thin-shell or linear objects (cloth, rope) with surface point clouds. Projective dynamics [23] were used to create a real-time physics-based model for tissue deformation, enhanced by a Kalman filter (KF) for refining the simulation with surface marker data. In a related work [24], they integrated FEM-based simulation into a deep reinforcement learning framework for grasping point policy learning, relying on offline stereo calibration for registration with the real world. For further improvements, [25] utilized a variational autoencoder with graph-neural networks to learn low-dimensional latent state variables' probability distributions. These variables were iteratively updated using an ensemble smoother with data assimilation to align the simulation with real data. However, their offline training for specific FEM simulation datasets poses challenges in realworld applications, especially in surgery, where lengthy datacollection and pretraining phases are impractical.
In this work, we propose an online simulation parameter optimization framework for bridging the real-to-sim gap in deformable object manipulation. We believe it offers significant benefits in recovering a realistic physics simulation model for planning and control in DOM scenarios, including surgical autonomy. To this end, we present the following novel contributions:
• An online "real-to-sim" residual mapping module is seamlessly incorporated into a physics simulation loop to predict a residual deformation ∆ t over time, which modifies simulation states to align with point cloud observations. • The module estimates geometry-aware deformation for sub-surface simulation particles by considering geometric constraints. Therefore, it can be applied to thin-shell and volumetric deformable objects, covering wider surgical use cases than previous works.
• An online simulation optimization framework is used that adaptively updates stiffness parameters through loss functions that are informed by the residual mapping module, enhancing the simulation's ability to predict future deformation and minimizing the the "real-to-sim" gap.
Let the state of a soft body in a physics simulation at time stamp t be x t ∈ R n×3 , where n is the number of particles of a mesh representing the soft body. We use the surface point cloud projected from the stereo-depth estimation as the observation, denoted as z t ∈ R m×3 , where m is the number of points in the point cloud observation. Let u t be a point-based positional control that is applied to the real object and simulation simultaneously. In this work, we use an extended position-based dynamic (PBD) simulator, which formulates constraints with positional and geometric data. A constraints-based formulation of PBD is given as
where Γ are static (boundary) conditions that are enforced at each simulation step. The set of geometric constraints that define the deformation across the simulation is
is the set of weighting parameters associated with any kind of non-boundary constraint. One can interpret that {C, k c } jointly defines the total energy potential (unitless) generated by constraints in a simulation:
This is used by PBD to iteratively update particle positions, by minimizing this energy term:
in which M is a diagonal mass matrix and λ is a Lagrange multiplier vector [26]. Note, while we demonstrate the proposed method with PBD simulation, it can be extended to real-to-sim problems in other differentiable simulators (e.g., FEM, projective dynamics [27]) without losing generality. We are only including the simulation step for PBD so one can follow the use of parameters through a forward simulation.
In this work, we aim at optimizing stiffness parameters associated with each particle, k ∈ R n . This converts to elastic constraints' weights in the simulation, denoted by k d ¢ k c , by averaging stiffnesses across all involved particles. The value of k d,i ∈ k d , that is, the weights of an elastic constraint
where particle indices Q i are considered by C i (x). Here, Q i contains card(Q i ) (cardinality) number of particles that Fig. 2. A flow chart of the proposed residual mapping module in the simulation loop. At each time step, a control u t is applied to both the real tissue and the simulation. In response, the PBD simulation solves for x t . A perception pipeline processes imagery data to obtain a surface point cloud z t of the tissue. The residual mapping module estimates the residual deformation ∆t via optimization, which is then used to update the simulation state.
are connected via mesh edges, triangles, and tetrahedrons. This achieves non-homogeneous elastic stiffness across different regions of a soft body. Given that we generate a discretized mesh solely at the initial step, without engaging in any re-meshing procedures throughout the simulation, this imposes challenges on accurately representing isotropic material properties within the tissue. In light of this, we propose considering k as a spatially-variable stiffness.
We apply the developed framework on two types of deformable bodies: thin-shell objects that are simulated with a single layer of particles (similar to cloth) and volumetric objects that are simulated with tetrahedral meshes.
We will learn a residual mapping module that is incorporated into the simulation loop to estimate the real-to-sim gap. This module will enable matching between simulation particles to a point cloud observation as shown in Fig. 2. The method employs a gradient-based non-rigid point cloud registration to estimate the residual deformation ∆ t ∈ R 3×n . To make the predicted residual mapping accurate and physically realistic, we consider both point cloud similarity D(•) and physical realness E(•) as cost functions to minimize. Both metrics are discussed in depth later.
Because the correspondence between the simulation particles and point cloud observation is unknown, we use Chamfer Distance as a measurement of similarities between two point clouds, which is computed by summing the squared distances between the nearest neighbor of two point clouds. For a thin-shell simulation, it is defined as
A volumetric mesh can't be directly aligned to a surface point cloud observation by minimizing Equation 4 as internal particles are not observed. Therefore, we alternatively minimize D(x ′ t , z t ) where x ′ t ∈ x t is the surface particles of the volumetric mesh, which is available from mesh initialization.
We enforce physical realness E(•) utilizing the simulator's geometric constraints. This is achieved by directly minimizing the simulation's energy potential:
where k ′ c is an uniform user-defined weight matrix and C are PBD constraints. Note that because sub-surface particles are not observable for volumetric meshes, this term is also essential for inferring their deformation. A visualization of how D(•) and E(•) affect the mesh deformation is shown in Fig. 3. Finally, a residual mapping function g(x t , z t ) is defined as
The above minimization problem is solved by performing 30 gradient descent steps with a learning rate of 50. At the end of a simulation step, ∆ t is used to update the current simulation state x t . Later, ∆ t is also viewed as a proxy to the real-to-sim gap that we seek to minimize through an online stiffness optimization approach.
The proposed online optimization method differs from previous real-to-sim methods as it does not rely on training on previously collected trajectories. In this way, we are solving an online problem that is much more generalizable to real-world, unstructured environments. The algorithm is also embedded in the simulation loop as shown in Alg. 1.
One term we want to minimize directly is the residual gap, which characterizes how much the simulation deviates from the observation for the current time step.
Its partial derivative with respect to stiffness parameters k is
where the first term is available by differentiating through the residual mapping module, and the second term is obtained by differentiating through a simulation step in Equation 1.
The residual gap alone doesn't consider historical information and, therefore, is sensitive to current observation noise. // Initialize a mesh and constraints.
To address that, we introduce a history term computed over a set of previous time points H t . It is defined as
where x h and ∆ h are snapshots of simulation states and residual mapping at time point h. Specifically, H t is constructed by uniformly sampling four snapshots from a window of the closest 20 previous frames. This loss function finds the stiffness parameters that keep each snapshot at rest when no control is provided (i.e., balancing external forces and internal elastic forces). This is a reasonable assumption when the system is not moving quickly.
In addition, we encourage a smooth spatial stiffness distribution by penalizing k's differences between neighboring particles. Let F be a set of all faces or tetrahedrons that sorts tuples of particle indices, the smoothness loss is written as
where k i ∈ k is the stiffness value of the i-th simulation particle. Note that intersection of different types of tissues is more complicated. which may be addressed by encouraging smoothness within individual semantic contours [28]. Lastly, all terms are summed up and back-propagated to the stiffness parameters. They are updated by taking a stochastic gradient descent step in every simulation step.
Fig. 4 shows our physical experiment setup. Our experimental procedures involve the utilization of the da Vinci Research Kit (dVRK) [29], employing its robotic gripper to precisely manipulate soft tissue by executing predefined trajectories. The trajectories are shown in Fig. 4 (e). Simultaneously, a stereo reconstruction pipeline processes stereo images captured by the da Vinci endoscopic camera in 720p, producing tissue surface point clouds. A piece of chicken skin and chicken muscle are used as subjects in our thin-shell and volumetric experiments, respectively. They are roughly 3mm and 2 cm thick. In both experiments, we used metal pins to fix one side of the tissue onto the bottom plane.
We utilize the Raft-Stereo [30] for stereo disparity estimation. Segment-Anything [31] aids in identifying image pixels corresponding to the tissue, allowing us to extract the tissue's surface point cloud from depth images. Subsequently, we employ inverse camera projection to convert the segmented depth into 3D positions. We also employ ArUco markers to determine the camera-to-world transformation. The point cloud is down-sampled to a size of 9000 points. Meshes are reconstructed by applying a ball-pivoting algorithm on the initial surface point cloud observation. For volumetric experiments, we process the surface meshes by adding thickness using the solidify modifier in Blender. Meshes are remeshed to contain 600 particles. The simulation's boundary conditions are selected at the locations of pins. The proposed residual mapping module requires a knowledge of surface particles that are visible from the camera. While we manually selected them, it is possible to select them automatically through visibility culling. The robotic gripper trajectories are manually labeled on the image.
We collect in total four trajectories, two for thin-shell and two for volumetric experiments. Later, they are referred as Thin-shell-1, Thin-shell-2, Volumetric-1, and Volumetric-2.
This framework is implemented in Pytorch [32]. Our PBD simulation adopts three types of geometric constraints: distance, volumetric, and shape-matching. They respectively preserve the distance between connected simulation particles, the volume of tetrahedral, and shapes formed by neighboring particles. We use distance and shape-matching constraints when simulating thin-shell tissue, whereas all three are used for volumetric tissue. For consistency, we optimize distance With the residual mapping, the simulation mesh matches to textured point cloud observation better than the original PBD. and shape-matching parameters for both thin-shell and volumetric test cases while keep volumetric constraint parameter fixed as k vol = 1e10. In addition, we impose bounds on our simulation parameters with a sigmoid operation. These bounds are k dist ∈ [0, 10] and k shape ∈ [0, 0.02], which are empirically found effective to optimize in. As optimizationbased methods are known to be sensitive to initialization, we evaluate the proposed method with three sets of initial parameters. They are k 1 : {k dist = 5, k shape = 0.15}, k 2 : {k dist = 1, k shape = 0.1}, k 3 : {k dist = 0.2, k shape = 0.005}. All stiffness parameters are uniformly initialized. Note that initializations can be determine from prior knowledge of the deformable object. Here, we aim at testing every methods thoroughly. They are optimized with an Adam optimizer with a learning rate of 0.1. Later, we refer to a PBD simulation with a residual mapping module as PBD-RM and one that further performs online updates as PBD-RM-ON.
In this section, we evaluate the effectiveness of applying the proposed residual mapping. The simulations are initialized with stiffness k C 2 . We measure the error between simulation and observation with Chamfer distances. Fig. 5 shows that the errors are notably reduced in all trajectories using the proposed residual mapping module. They are stably kept at low levels throughout all trajectories. In the thin-shell cases, the proposed module shrinks the average errors from 1.84 mm and 2.55 mm to 0.001 mm and 1.31 mm. Likewise, the average errors are reduced from 4.24 mm and 5.79 mm to 1.78 mm and 1.93 mm in the volumetric case. Fig. 6 visualizes the reduced differences between point clouds and mesh with our mapping module. For the volumetric case, the module simultaneously aligns the volumetric mesh's surface particles to an observation and also deforms sub-surface particles, respecting PBD's geometric constraints.
Here, we evaluate the method's future deformation prediction capabilities. Specifically, we formulate an average future gap e t and a future keypoint error f t as our metrics. Knowing the current simulation state x t and a future control sequence, we can roll out a simulation, without knowing future observations, to get future state sequence [x t+1 ...x t+T ]. The average future gap is computed as
Let p t be a vector of keypoints' positions at time t, and let x i t represents the i-th particle in x t . The future keypoint error is defined as
TABLE II. COMPARISON OF TIME-AVERAGE (et, ft) OF DIFFERENT METHODS WITH PRE-DEFINED INITIALIZATIONS. RESIDUAL MAPPING (PBD-RM) AND ONLINE OPTIMIZATION (PBD-RM-ON) IMPROVE PREDICTION CAPABILITIES. Methods k 1 k 2 k 3 k 1 k 2 k 3 Thin-shell-1 Thin-shell-2 PBD 0.52, 2.24 0.47, 2.14 0.34, 1.74 1.06, 3.80 0.96, 3.63 0.89, 3.51 PBD-DiffCloud 0.45, 1.76 0.46, 1.73 0.33, 1.57 0.90, 4.22 0.97, 4.14 1.02, 4.28 PBD-RM 0.46, 2.23 0.42, 2.13 0.32, 1.77 0.84, 3.65 0.76, 3.47 0.67, 3.23 PBD-RM-ON 0.40, 2.04 0.31, 1.66 0.30, 1.72 0.65, 2.98 0.58, 2.72 0.62, 2.89 Volumetric-1 Volumetric-2 PBD 1.20, 3.15 1.59, 4.03 2.51, 6.57 1.75, 4.45 2.46, 5.53 3.33, 7.58 PBD-RM 1.02, 3.06 1.13, 3.52 1.30, 4.32 1.41, 4.33 1.65, 5.00 1.86, 5.70 PBD-RM-ON 0.94, 2.91 0.97, 3.21 1.00, 3.38 1.40, 4.22 1.48, 4.56 1.52, 4.71 Here, keypoint displacement ∆p is approximated by the displacement of p's closest neighboring particles at t = 0.
For all experiments, we pick a future horizon T = 10. A set of 15 keypoints are manually labeled on images every 10 frames for each trajectory. These keypoints are in deformed regions and are selected based on their visual features.
For thin-shell experiments, we compare to DiffCloud [21], a baseline that optimizes simulation parameters by minimizing point cloud differences on a training trajectory. we optimized with the first 30 frames of each trajectory for 25 epochs. Table II shows a comparison on the two aforementioned metrics. Compared to the original PBD, PBD-RM shows reduced errors in all experiments, indicating that correcting the current simulation state with residual deformation leads to better future prediction. Compared to DiffCloud, PBD-RM-ON doesn't consistently outperform it on the relatively shorter and simpler Thin-shell-1. However, it performs significantly better than all other methods in the other three experiments. DiffCloud produces unsatisfactory results on Thin-shell-2 as it overfits to the beginning of that trajectory, whereas the proposed method does not as it updates parameters online. Fig. 7 supports the same observation, visualizing the average future gaps over time. Our PBD-RM-ON largely reduces the gap in the middle of trajectories where deformations are Fig. 8 visualizes the spatial distribution of future gaps. It indicates larger errors located where tissues are experiencing bending. Both PBD-RM and PBD-RM-ON are effective at reducing errors in those regions. In Fig. 9, the spatial distribution of optimized parameters is overlaid on real images. The stiffness parameters exhibit spatial and temporal variations. Effects such as distance stiffness lowering around the deforming areas in the thin-shell case can be attributed to the chicken skin's inherent softness. In cases of heterogeneous parameter distributions, we hypothesize that their spatial variation may better approximate the complex behavior of the real tissue. However, whether the spatial variations represent can real material stiffness needs further future studies.
In this work, we have introduced an online simulation optimization framework designed to address the real-to-sim gap that arises in deformable tissue manipulation. A residual mapping module is seamlessly integrated into a simulation loop, achieving a minimal Chamfer distance between simulated particles and observation while preserving the geometric relationships in the simulator. Our optimization approach updates constraints' stiffness parameters online. In real tissue experiments, it is proven to be effective at improving predictive performance. This improvement holds the potential for enhancing model-based control for DOM, such as identifying an accurate simulation model for model predictive control in autonomous surgeries. For deploying this work into real applications, one limitation is its computation speed. Currently, PBD-RM-ON takes 0.9s and 2.6s to complete a step for thin-shell and volumetric meshes. Additional, our perception pipeline takes 0.4s per image. Despite that, it can be accelerated by GPU implementation or asynchronous computation. Other key issues, such as perception failure due to noise and unmodeled phyiscal effects such as other external foces, need to addressed in future works. A possible future avenue of this work would be learning residual non-linear constraints to capture more intricate tissue behaviors.