Efficient Penetration Depth Computation between Rigid Models using Contact Space Propagation Sampling
Abstract
We present a novel method to compute the approximate global penetration depth (PD) between two nonconvex geometric models. Our approach consists of two phases: offline precomputation and runtime queries. In the first phase, our formulation uses a novel sampling algorithm to precompute an approximation of the highdimensional contact space between the pair of models. As compared with prior random sampling algorithms for contact space approximation, our propagation sampling considerably speeds up the precomputation and yields a high quality approximation. At runtime, we perform a nearestneighbor query and local projection to efficiently compute the translational or generalized PD. We demonstrate the performance of our approach on complex 3D benchmarks with tens or hundreds of thousands of triangles, and we observe significant improvement over previous methods in terms of accuracy, with a modest improvement in the runtime performance.
I Introduction
Accurate and efficient computation of interpenetration is important in many areas, including computer graphics, haptics, and robotics. A common metric that is used to measure the extent of interpenetration between two intersecting objects is penetration depth (PD), which is defined as the minimum amount of movement or transformation required to separate two incollision objects. The resulting motion may correspond to translational alone (translational PD) or to both translational and rotational motion PD (generalized PD). PD computation is frequently used for many applications, such as physicallybased simulation [1], samplebased motion planning [2], haptics [3, 4], and contact manipulation [5].
Computing the exact PD in 3D is a challenging task because of the time complexity involved in translational PD and the worst case time complexity for generalized PD, where and are the number of triangles in two nonconvex input models [2]. Given the high combinatorial complexity of exact PD computation, many approximate algorithms have been proposed. Some of the simplest algorithms compute the intersecting features of these two models and use them to compute local PD that is based on a measure of separating those overlapping features. In fact, current game engines such Box2D [6] and Bullet [7] use local PD computations for collision response. However, the accuracy of local PD algorithms depends on relative configuration of two objects [8, 9]. Other techniques are based on computing an approximation of the configuration space boundary [10, 11, 12], but the accuracy of these techniques can vary for different configurations of two objects and it is hard to derive tight bounds. There are no good reliable algorithms for global PD computation between arbitrary nonconvex 3D shapes.
Main Results: In this paper, we present a novel algorithm to approximate global PDs between rigid objects based on efficient sampling in the contact space. Our approach can compute both translational and generalized PD with high accuracy for nonconvex models. We first precompute an approximation of the contact space of two overlapping objects by generating samples in the contact space. We generate our initial samples using random sampling and use a novel propagation algorithm to generate additional samples via local search. The use of propagation sampling considerably speeds up precomputation and results in a high quality contact space approximation. At runtime, our algorithm performs a nearestneighbor query to compute the PD. We also analyze the properties of our sampling scheme and highlight its benefits. Compared with prior PD algorithms, our approach offers the following benefits:

The overall algorithm is general and directly applicable to complex nonconvex and nonmanifold models.It can compute translational and generalized PDs.

The use of propagation sampling can considerably accelerate the precomputation and provides a high quality approximation of the contact space.

The runtime query is very fast (a few milliseconds) and can be used for interactive applications.

The overall algorithm is more accurate as compared to prior local and global PD computation algorithms.
We highlight the performance of our algorithm on different models, which contain tens or hundreds of thousands of triangles with sharp features. We also highlight the considerable improvements in the accuracy of the runtime query compared with recent algorithms based on active learning [10] and local optimization [13, 4]. In particular, our approach can considerably reduce the error in PD computations over these methods. This paper is an extension of our work [14] and provides additional technical details that were not included in [14] due to the space limitations.
The remainder of this paper is organized as follows. In Section II, we survey the literature related to the configuration space and PD computation. We introduce our notation and give an overview of the algorithm in Section III. We present our contact space sampling algorithm in Section IV, and we discuss and analyze its performance on many complex benchmarks in Section V and Section VI.
Ii Related Work
Iia Configuration Space Computation
There is extensive work on configuration space computations in robotics, geometric computing, and related areas. In the most general cases, configuration space computations can be reduced to computing the arrangement of contact surfaces [15]. However, these approaches are susceptible to robustness issues. Moreover, the worstcase complexity of the arrangement computation can be as high as , where is the number of contact surfaces in the arrangement and is the dimension of the configuration space [16]. Some techniques for approximating the configuration space in lower dimensions are based on generating a discrete number of slices [17]. When the object movement is limited to translational motion, the resulting configuration space corresponds to the Minkowski sum of two objects [18, 19].
A substantial amount of work in motion planning involves approximating the configuration space with sampling techniques. These include various randomized algorithms that compute roadmaps for collisionfree path planning. Some of these approaches, such as [20, 21, 22, 23], also consider the problem of sampling in a constrained configuration space or a manifold, which is similar to the contact space sampling discussed in this paper. However, these methods are designed for collisionfree motion planning, and hence only require the generated samples to capture the connectivity of the free part of the constrained configuration space.
IiB PD Computation
Given two convex polytopes, we can compute the exact translational PD using Minkowski sum computation [24, 25, 26]. For nonconvex polyhedral models, the PD can be computed using a combination of convex decomposition, pairwise Minkowski sums, and exact union computation [11]. Different techniques have been proposed to approximate the boundary of the Minkowski sum [11, 15], but they are limited to offline and noninteractive applications. Most practical algorithms for translational PD are based on local computations. These local algorithms only consider the intersecting or overlapping of features such as the vertices, edges, and faces. Based on pairwise intersections, they tend to estimate a motion that would separate these intersecting features [27, 9, 12, 28, 29, 30]. Other techniques estimate the local intersection volume and its derivative to perform volumebased repulsion [31]. Local translational PD computation can also be estimated using distance fields [8]. Pointbased Minkowski sum approximation [32] has been used approximate the translational PD. The exact computation of generalized PD can be formulated in terms of computing the arrangement of contact surfaces [2]. However, no practical algorithms are known for exact computation due to its high combinatorial complexity. Most practical algorithms are based on local optimization techniques [33, 34, 4, 13]. However, due to the high time and storage complexity, most generalized PD algorithms are based on local optimizationbased techniques [33, 34, 4, 13]. [10] recently proposed a learningbased approximate penetration depth computation algorithm that reduces the contact space problem to robust classification by finding a separating surface between incollision and collisionfree samples in the configuration space. However, this algorithm cannot provide high quality approximation of the contact space of objects with sharp features, because it represents the contact space using the SVMs (support vector machines). Recently, Kim et al. [35] present a hybrid PD computation algorithm that combines this active learning approach with local optimization based methods to improve its accuracy.
Iii Background and Overview
In this section, we introduce our notation and give an overview.
Iiia Contact Space
We denote the configuration space for a pair of triangular meshes and as . Each configuration or point in the configuration space represents the relative transform (i.e., position and orientation) of with respect to . In the rest of this paper, we assume that is movable and is fixed. The configuration space is composed of two parts: collisionfree space represented as , and incollision or obstacle space represented as , where corresponds to located at the configuration , and and correspond to set closure and interior operations, respectively.
The boundary of is called the contact space and is denoted as . The contact space corresponds to the configurations where and just touch each other without any penetration. Moreover, a contact configuration is classified as a collisionfree configuration in our formulation.
IiiB PD Formulation
The global penetration depth corresponds to the minimum motion or transformation required to separate two intersecting objects and [25, 26]:
(1) 
where corresponds to an incollision configuration and is a configuration that belongs to the contact space . We use the notation to represent a distance metric between two configurations. This includes the Euclidean metric for translational PD, and many different formulations can be used for generalized PD computation. We denote as the contact configuration where PD achieves its minimal value: .
Different formulations of PD can be defined by appropriate metrics. The metric for the translational motion () is simple and is the standard Euclidean distance metric between vectors corresponding to the configurations. The metric for the general motion () can be defined using different formulations, including the weighted Euclidean distance [3], object norm [36, 4, 13], and a displacement distance metric [34]. In this paper, we use the object norm [36, 4, 13] as the metric, which can be intuitively defined as an average squared length of all displacement vectors between two objects.
IiiC Approximate Computation
In order to construct an approximation of the contact space, we perform an offline sampling in the configuration space, as shown in Figure 1. Given two input objects and , our method starts with a randomsample in the contact space. This randomsample can be generated using traditional continuous collision checking (CCD) techniques. These techniques compute the first time of collision or contact between two objects by reducing the problem to finding roots of polynomial equations corresponding to the triangle features. Continuous collision detection has been widely used for physicallybased simulation [29] as well as local planning in robotics. Next, our method performs an iterative local search around the initial randomsample by sliding object over object ’s surface, and generates more samples on the contact space. We denote the samples generated on the contact space during the local search as propagatesamples to distinguish them from randomsamples (see Figure 1(h)). The local search stops when no more propagatesamples can be generated, and we then restart the iteration with a new randomsample in the contact space. This iterative process continues until a sufficient number of samples have been generated. The randomsamples and propagatesamples computed by our approach make up an approximate samplebased representation of the contact space between and , and we denote this approximation as .
IiiD Approximate PD Computation
Given the approximate representation of the contact space, , we compute the approximate global PD by performing a nearestneighbor query in . The definition of approximate penetration depth is analogous to the exact penetration depth in Equation 1:
(2) 
where the domain for is restricted to . The accuracy of is governed by the accuracy of with respect to . Given a query configuration , we perform a nearestneighbor search to find the configuration that is closest to the decision boundary . Finally, the distance between and is computed using an appropriate distance metric and the result is an approximation of the exact PD value. As mentioned in Section IIIB, we use the object norm as the distance metric.
Iv Contact Space Propagation Sampling
In this section, we present our contact space propagation sampling algorithm that computes an offline approximation of the contact space. Our sampling algorithm is an iterative algorithm. During each iteration, we start from a randomsample in the contact space, and then perform a local search around this initial sample to generate more propagatesamples on the contact space. Once the local search stops, we repeat the iterative step with a new randomsample. This iterative process continues until a sufficient number of samples have been generated. The randomsample on the contact space during each iteration is computed by first generating two samples in the configuration space, one in collisionfree space and the other in the incollision space. We join those samples by a straight line in the configuration space and find its intersection with the contact space. This reduces to computing the first time of contact between a collisionfree and an incollision configuration, which corresponds to a CCD query. The resulting sample on the contact space is the randomsample used during this iteration.
Ideally, the local search procedure should run many steps and generate sufficient numbers of propagatesamples, in order to cover a large portion of the contact space. This coverage is important for the efficiency of the sampling algorithm, because the generation of a randomsample requires the expensive CCD query. This query is more expensive than the generation of a propagatesample that only needs to perform the DCD (discrete collision detection) query. If the local search can generate a high number of propagatesamples, it amortizes the computational cost of generating a randomsample over a high number of propagatesamples, and improves the efficiency of our precomputation step. Our main goal is to design a fast and effective local search algorithm that can compute the propagatesamples quickly. To that end, we perform a breadthfirst propagation on the contact space, starting from the randomsample. The breadthfirst propagation maintains a queue (we call it the propagatequeue) of contact samples. During each step of this propagation, we pop one contact sample from this queue, and then slide object over the surface of object in different directions along its boundary. This propagation step results in a set of new samples around , as shown in Figure 1. Next, we perform collision checking for these samples , and only add the collisionfree samples into the propagatequeue. The breadthfirst search is repeated until the queue is empty.
A new sample is propagated from a contact sample and is added into the propagatequeue only if it is collisionfree. Otherwise, the sample would be discarded and the actual execution of the local search’s propagation may be interrupted, as shown in Figure 2. To overcome this challenge, we classify the local search process into two cases: the boundary configuration case when , and the internal configuration case when . In the internal configuration case, we resume the local search computation according to objects’ contact features (i.e., local vertices, edges, faces). This formulation greatly improves the efficiency of the local search. The overall local search computation algorithm is shown in Algorithm 1.
Iva Boundary Configuration Case
In the boundary configuration case, each step of the local search is a standard propagation step from a contact sample . We first compute the contact pair between and , where and are two contact points on objects and . The points and belong to two different objects but overlap with each other, just touching at the configuration . We also compute the angle between the normal vectors at and . Next, we slide the object over the surface of object with a distance , which can be any value less than the edge length of triangular mesh(since we generate the samples on the mesh of object). For our case, we we use the edge lengths of the fixed objects as the sliding step, and thus the propagatesamples all land in the vertices of object . During the sliding movement, the contact point on remains unchanged as , and the contact point on moves from to (see Figure 2). Now and touch at the new contact point pair . We further rotate the object such that the angle between the two objects’ contact normals remains to be (as shown in Figure 2). In this way, we compute a new configuration , which can be specified using , and for 2D objects. Similar propagation procedure can also be defined for 3D objects, whose configuration space has 6 dimensions.
We represent the sliding movement from to as a transition function . The new generated configuration is pushed into the propagatequeue for future propagationsample computations. This sliding procedure is executed along different directions on the surface of object around , and results in a set of new configurations spreading over the neighborhood of on the contact space. The collisionfree samples in are located on the contact space and we add them directly into the propagatequeue. For incollision samples in , we treat them as the internal configuration case and stop propagation.
In the above description, we restrict all new propagatesamples to have the same angle between the two objects’ contact normals. This simple heuristic greatly increases the probability that a new sample will be collisionfree. It is based on the assumption that the surface curvature around the neighborhood of is roughly constant. Therefore, a configuration with the same relative angle as should have a high probability to be collisionfree after the sliding movement.
The parameter in a boundary configuration propagation step determines the stepsize of the propagation. Its value varies during the local search process. In particular, is inversely proportional to the relative scale of with respect to , and it is also related to the surface curvature at .
IvB Internal Configuration Case
In this case, is inside , i.e. inside the Cobstacle space. This can happen when two objects are very close in size or the surface of is ‘bumpy’, i.e., the curvature changes dramatically over the surface. In this case, the propagation search step usually explores only a few steps because and will collide even when only slides a small step over ’s surface. In an extreme situation, every local search returns no propagationsamples and all samples generated on the contact surface are randomsamples. This will result in very slow sampling procedure, and these samples cannot be evenly distributed over the contact space.
Our solution for the internal configuration case exploits the contact features, and is based on the following property of the sliding movement:
Theorem 1
Suppose one step of slide moves a contact sample into an incollision sample . The transition function is
(3) 
where is the contact pair at , is the contact pair at , and is the angle between contact boundary configurations. On the resulting sliding trajectory, there exists one configuration such that and are in contact, but have at least one additional contact point other than . Here we assume the sliding movement is parameterized by . We call the configuration the critical configuration.
For each , we denote as its position corresponding to the configuration . We further define , i.e., the first time that contacts with on a point other than . Since is collisionfree and is incollision, we know . has the property that contacts with , but there is at least one additional contact point other than . Therefore, the theorem is proved.
Based on this theorem, we resume the propagation at an incollision sample by first computing the critical configuration between and . has more than one contact point, and we denote the set of contact points other than as . For each contact point in , we continue the local search step by changing the contact point on that remains unchanged during the movement from to . In particular, we compute the contact pair and , the angle between the contact boundary configurations at and . Object now starts sliding from the contact point and the corresponding transition function is . Figure 3 illustrates this process.
To avoid generating repeated samples during the internal configuration case, we use a kdtree to check the new contact configuration and use it as the new propagation seed only if it is not close to any existing samples, as shown in the line in Algorithm 1.
IvC Runtime PD Queries
In the PD query stage, we use the contact configurations generated during the precomputation. Given a query configuration , we use nearestneighbor query to find two contact configurations and that are closest to and incident to a same triangle on the fixed object. Next, we compute a configuration which is a linear combination of these two contact configurations, and the line connecting an is perpendicular to the line connecting and :
(4) 
We then perform a linear search from along directions of until we find a collisionfree configuration . The nearestneighbor query is performed based on the metrics listed in Section III. Finally, we have to compare the distance of with those of , , and choose the smallest one as the PD value of .
V Performance and Comparison
In this section, we highlight the performance of our propagation sampling based precomputation and the runtime PD query on a set of challenging 3D benchmarks as shown in Figure 4. We use PQP query package to perform all the collision tests used in our framework. For PD query, we investigate both the translational and generalized PD, using the distance metric as mentioned in Section III.
We implement our algorithm in C++ on an Intel Core i7 CPU running at 3.30GHz with 16GB of RAM on Window 7 (64bits) PC. All the performance and timing results are generated using a single core. Our precomputation algorithm can be easily parallelized on a multicore PC because both the randomsamples and the propagatesamples on the contact space can be generated in parallel.
Va Performance
Time Cost: We now highlight the performance of both the offline precomputation and runtime PD queries. Table I and Table II show the time costs of the precomputation phase for general PD and translational PD respectively, given different number of propagation steps. We use at most 100,000 propagation steps, and each step will generate one randomsample and a set of propagatesamples extended from the randomsample. We can see that the time cost for a single propagation step increases while more samples are generated, and this is due to the fact that while more samples are generated, it is more difficult to find unrepeated new samples in the contact space. Table I and Table II also provide the number of samples generated for each propagation step, which is averaged over all the 100,000 propagation steps. Table III provides our method’s runtime cost for a single PD query, which is computed as the average time cost after computing 1,000 randomly generated incollision queries.
#propagations  30,000  50,000  80,000  100,000  #samples 

Donut  314  551  693  891  79 
CAD1  421  569  745  883  411 
CAD2  390  581  737  923  323 
Teeth  520  892  1,455  1,621  1ï¼231 
Dragon  590  940  1,401  1,632  2ï¼997 
Buddha  603  875  1,337  1,501  1ï¼820 
#propagations  30,000  50,000  80,000  100,000  #samples 

Donut  156  193  344  423  63 
CAD1  159  243  387  659  329 
CAD2  163  247  401  613  341 
Teeth  155  471  912  973  768 
Dragon  333  752  941  1,429  667 
Buddha  346  771  1,017  1,918  751 
Convergence: We investigate the convergence of our sampling algorithm by varying the number of samples and evaluating its benefits in terms of approximating the contact space. Since it is difficult to obtain the ground truth for the general PD, the investigation is only for the translational PD. We use the result from the Minkowski sum based approach [32] as the ground truth of the contact space, and then we compute the error between the PD results from the ground truth and the PD results provided by our propagation sampling approach.
As we generate more randomsamples and propagatesamples, our precomputation algorithm provides a better approximation of the contact space. To evaluate the approximation quality of the contact space, we measure the number of vertices in the original models that have been visited during the propagation as contact points (e.g., and in Figure 2). This is due to the fact that the denser the contact space is sampled, the more vertices should be visited as contact points during the sampling process. The measurement result is shown in Figure 5, and we can observe that for all benchmarks, the propagation procedure convergence after generating million samples, and a high quality approximation to the contact space is achieved.
VB Comparison with Prior Methods
We compare the performance of our algorithm with two recent methods that can compute translational and generalized PD. One is a global approach [10] that uses active learning to generate the samples in the entire configuration space, and then computes the contact space approximation as the surface separating collisionfree and incollision samples. Its approximation of the contact surface is relatively smooth due to its leverage of SVMs and hence this method may not provide accurate PD results for query configurations with small penetrations. Overall, this approach provide a conservative bound on the PD. The second method uses local optimization [13, 4] to compute translational and generalized PD. It starts from a good initial guess to the PD result and optimizes the computation based on appropriate metric. Its accuracy depends on the initial guess and it may get stuck in a local minimum. Unlike our method or the active learning approach, this method has no precomputation. Therefore this approach can also be used for PD computation involving deformable models. Recently, Kim et al. [35] presented a hybrid approach that combines the global method [10] with the local method [13] to overcome some of these problems.
Runtime Performance: We compare the performance of runtime query of our approach with these methods, and the result is shown in Table III, for both general and translational PD formulations. We observe that our approach is significantly faster than other two approaches, because our method can generate samples more evenly and densely distributed in the contact space, and thus query configuration can easily find a nearestneighbor contact configuration in fewer iterations.
Accuracy of PD computation: In our experiments, our method quickly generates more than 1 million samples in 1012 minutes with less than 5 MB storage. The error of our online PD query is about 5% of the actual PD value for both the translational and general PD, while this value is 1520% for [10] and 10% for [13, 4]. The main reason for this improvement is because our method searches over the entire contact space globally while the local optimization methods [13, 4] search along the gradient direction. The gradient search strategy can only generate suboptimal results, especially for complex models. Our approach also outperforms [10, 35] which approximates the contact space by active learning, because it can provide a more accurate approximation of the contact space.
We compare the accuracy between our method and the active learning approach and the local optimization approach by computing the perframe error on a set of different benchmarks (CAD1, Bunny, Lion, Dragon and Donuts). The result is shown in Figure 6. From the results, we can observe that the PD error of our method is about 5% for both the translational and general PD queries. The error is mainly caused by the resolution of the samples. But as we discussed in the algorithm, all the samples in our method is located on the contact space, and thus the result is much more accurate than that of [10] which fits the contact space by machine learning method. While both methods use the same number of samples for general PD computation, the error of [10] is about 10%  20%, as shown in the third row of Figure 6, while our method is about 5% as shown in the fourth row of Figure 6. For the translational PD, the error of our approach is about 3%, which is smaller the the corresponding results in local optimization approaches [13, 4] as shown in the second row of Figure 6. Local optimization approaches also search over the objects’ surface, but they do not perform well for objects with complex shapes. This is because they only search along the gradient direction and may only obtain suboptimal solutions. For instance, they may miss the optimal results in areas of sharp features. Our method will not miss these areas in most situations since we can generate 1 million samples over the contact space, which help our method to explore the entire contact space. By searching through these samples, our method can generate more accurate results than gradient searching method. Overall, our approach is more accurate than prior PD methods.
Donut  CAD1  CAD2  Dragon  Teeth  Buddha  
Our ()  0.93  0.95  1.03  1.18  0.98  1.31 
Our()  0.81  0.85  1.13  1.06  1.01  1.11 
Active [10]  0.95  1.41  1.32  1.43  0.92  1.48 
Local [13]  6.56  78.3  54.4  132  111  153 
Vi Analysis
In this section, we evaluate some properties of the propagation sampling algorithm used in our precomputation phase. In particular, we discuss 1) the bounds on the time complexity of our precomputation scheme, and 2) why our approach can generate samples with better distribution in the contact space than pure random sampling.
Via Time Complexity
For the time complexity analysis, we assume the object can only perform translational movements and has vertices, and object is a connected mesh with vertices, where . We also denote and as the time costs for one continuous collision checking and one discrete collision checking respectively. On average, in our 10 thousand random tests. Then we have the following estimation for the precomputation’s time complexity:
Theorem 2
The precomputation’s time complexity has a lower bound of , and has an upper bound of .
To obtain a high quality samplebased representation for the contact space, one way is to generate all the configurations where contacts with at one of ’s vertices. Ideally, we only need to generate one randomsample, and all the other configurations can be visited in one iteration of propagation process. The time complexity for such ideal case is , which is a lower bound of the propagation sampling’s time complexity.
The upper bound of the propagation sampling can be achieved by considering the expected time cost of a pure random sampling required to visit all the vertices of object . Suppose the sampling process has already visited vertices of , the expected time for the sampling process to find a new vertex different from the visited vertices is . As a result, the expected time cost for the entire sampling procedure is , which can serve as a upper bound of our propagation sampling’s time complexity.
The two bounds estimated in this theorem are conservative, but intuitively describes the performance of our precomputation scheme. In practice, the running time of our algorithm is much lower than the upper bound estimated, because for each propagation iteration started from a randomsample, we can generate a large number of propagatesamples due to our special treatment to the internal configuration cases in Section IV.
ViB Sample Distribution
We next show that our method can generate samples evenly distributed over the contact space. First, the propagation movement has the following property:
Theorem 3
Commutativity: Suppose the object slides along the surface of the object according to the transition function , where the transition function is as described in Equation 3. and are object ’s configurations before and after the transition, and and are the contact points on the object for these two configurations. Then the inverse slide can be formulated by the transition function .
This property is due to the fact according to our transition rule, the state and goal states of the transition function have onetoone correspondence. Intuitively speaking, this means that if the sample is visited first, then it can extend to using the transition ; if the sample is visited first, it can extend to using the same transition.
Next, we show that the samples generated by our propagation scheme will have no duplications:
Theorem 4
Uniqueness: Any two sets of samples generated by different propagation iterations will have no overlaps.
We use Figure 7 as an example for the two sets of sample configurations, where points belonging to different sets are marked as orange and blue respectively. Each sample set includes one seed randomsample (point and for this example), and a set of other samples extended from the seeds.
If the two sample sets indeed overlap with each other, and we assume the duplicated points are and , then according to Theorem 3, the point (also ) is able to propagate to , then and , and finally . This means that all the points in the second sample set will overlap with some points in the first set. However, while generating the randomsample for the second set, we use a Kdtree to guarantee that its distance to all points in the first set is large enough and thus such duplication for the randomsample is impossible. This theorem explains why our approach will not generate repeated samples and can result in an efficient precomputation phase.
Overall, our propagation sampling method has three benefits:

Efficiency Theorems 3 and 4 imply that no duplicated samples would be generated during the propagation. As a result, the algorithm will not waste time on generating redundant samples, and can make steady progress toward constructing a more precise representation of the contact space while generating more and more samples. These properties result in the error reduction shown in Figure 5. In addition, these two theorems also guarantee the even distribution of samples over the surface of the contact space.

Accuracy Our propagation steps directly slide the moving object over the surface of the fixed objects. As a result, every generated sample locates exactly on the contact space, and this results in the high accuracy of the PD results.
Vii Limitations, Conclusions and Future Work
We present a new PD approximation algorithm between general 3D models. We compute an approximation of contact space using propagation sampling. The propagation sampling scheme improves the accuracy of the approximation and is much faster as compared to only using random samples. We highlight the performance on many complex 3D models and highlight the benefits in terms of runtime performance and accuracy.
Our approach has some limitations. If the contact space has very narrow components, our sampling approach may miss them and thereby affect the accuracy of PD computations. It is possible that propagation sampling may not find samples during the local search computation. We would like to investigate the use of narrowpassage algorithms in samplebased motion planning to improve the performance. It would be useful to improve the performance of propagationsampling by utilizing the local curvature of the surface to choose appropriate directions. We would like to evaluate the performance on complex and articulated models, and also integrate with dynamics simulation and motion planning algorithms.
References
 [1] D. Baraff and A. Witkin, Physically Based Modeling. ACM SIGGRAPH Course Notes, 2001.
 [2] L. Zhang, Y. J. Kim, G. Varadhan, and D. Manocha, “Generalized penetration depth computation,” Computer Aided Design, vol. 39, no. 8, pp. 625–638, 2007.
 [3] D. Wang, S. Liu, X. Zhang, and J. Xiao, “Configurationbased optimization for six degreeoffreedom haptic rendering for fine manipulation,” IEEE Transactions on Haptics, vol. 5, no. 4, pp. 332–343, 2012.
 [4] C. Je, M. Tang, Y. Lee, M. Lee, and Y. J. Kim, “Polydepth: Realtime penetration depth computation using iterative contactspace projection,” ACM Transactions on Graphics, vol. 31, no. 1, pp. 5:1–5:14, Feb. 2012.
 [5] M. Koval, N. Pollard, and S. Srinivasa, “Manifold representations for state estimation in contact manipulation,” in International Symposium on Robotics Research. Springer, 2013.
 [6] E. Catto, “Box2D: A 2d physics engine for games,” http://box2d.org, 2010.
 [7] E. Coumans, “Bullet physics library,” http://bulletphysics.org, 2010.
 [8] B. Heidelberger, M. Teschner, R. Keiser, M. MÃ¼ller, and M. H. Gross, “Consistent penetration depth estimation for deformable collision response,” in International Fall Workshop on vision, modeling and visualization, 2004, pp. 339–346.
 [9] S. Redon and M. C. Lin, “A fast method for local penetration depth computation,” Graphical Tools, vol. 8, no. 1, pp. 63–70, 2006.
 [10] J. Pan, X. Zhang, and D. Manocha, “Efficient penetration depth approximation using active learning,” ACM Transactions on Graphics, vol. 32, no. 6, pp. 191:1–191:12, Nov. 2013.
 [11] Y. J. Kim, M. A. Otaduy, M. C. Lin, and D. Manocha, “Fast penetration depth computation for physicallybased animation,” in Proceedings of SIGGRAPH/Eurographics Symposium on Computer Animation, 2002, pp. 23–31.
 [12] J.M. Lien, “A simple method for computing Minkowski sum boundary in 3D using collision detection,” in Algorithmic Foundation of Robotics VIII, ser. Springer Tracts in Advanced Robotics. Springer Berlin / Heidelberg, 2009, vol. 57, pp. 401–415.
 [13] M. Tang and Y. J. Kim, “Interactive generalized penetration depth computation for rigid and articulated models using object norm,” ACM Transactions on Graphics, vol. 33, no. 1, pp. 1:1–1:15, Feb. 2014.
 [14] L. He, J. Pan, D. Li, and D. Manocha, “Efficient penetration depth computation between rigid models using contact space propagation sampling,” IEEE Robotics and Automation Letters, 2016, to appear.
 [15] G. Varadhan, Y. J. Kim, S. Krishnan, and D. Manocha, “Topology preserving approximation of free configuration space,” in Proceedings of International Conference on Robotics and Automation, 2006, pp. 3041–3048.
 [16] M. Sharir, “Algorithmic motion planning,” in Handbook of Discrete and Computational Geometry, J. E. Goodman and J. O’Rourke, Eds. Boca Raton, FL: CRC Press LLC, 2004, ch. 47, pp. 1037–1064.
 [17] E. Sacks and C. Bajaj, “Sliced configuration spaces for curved planar bodies,” International Journal of Robotics Research, vol. 17, pp. 639–651, 1997.
 [18] L. J. Guibas and R. Seidel, “Computing convolutions by reciprocal search,” Discrete & Computational Geometry, vol. 2, no. 1, pp. 175–193, 1987.
 [19] T. LozanoPérez, “Spatial planning: A configuration space approach,” IEEE Transactions on Computers, vol. C32, no. 2, pp. 108–120, 1983.
 [20] X. Ji and J. Xiao, “On random sampling in contact configuration space,” in Workshop on Algorithmic Foundations of Robotics, 2000.
 [21] L. He and J. V. D. Berg, “Efï¬cient exact collisionchecking of 3d rigid body motions using linear transformations and distance computations in workspace,” in IEEE International Conference on Robotics and Automation, 2014.
 [22] N. M. Amato, O. B. Bayazit, L. K. Dale, C. Jones, and D. Vallejo, “OBPRM: An obstaclebased prm for 3d workspaces,” in Workshop on the Algorithmic Foundations of Robotics, 1998, pp. 155–168.
 [23] O. Salzman, M. Hemmer, B. Raveh, and D. Halperin, “Motion planning via manifold samples,” Algorithmica, vol. 67, no. 4, pp. 547–565, 2013.
 [24] G. van den Bergen, “Proximity queries and penetration depth computation on 3D game objects,” in Game Developers Conference, 2001.
 [25] P. K. Agarwal, L. J. Guibas, S. HarPeled, A. Rabinovitch, and M. Sharir, “Computing the penetration depth of two convex polytopes in 3D,” in Proceedings of Scandinavian Workshop on Algorithm Theory, 2000, pp. 328–338.
 [26] Y. J. Kim, M. C. Lin, and D. Manocha, “DEEP: Dualspace expansion for estimating penetration depth between convex polytopes,” in Proceedings of International Conference on Robotics and Automation, 2002, pp. 921–926.
 [27] E. Guendelman, R. Bridson, and R. Fedkiw, “Nonconvex rigid bodies with stacking,” ACM Transactions on Graphics, vol. 22, no. 3, pp. 871–878, 2003.
 [28] M. Tang, M. Lee, and Y. J. Kim, “Interactive hausdorff distance computation for general polygonal models,” ACM Transactions on Graphics, vol. 28, no. 3, pp. 74:1–74:9, 2009.
 [29] M. Tang, D. Manocha, M. A. Otaduy, and R. Tong, “Continuous penalty forces,” ACM Transactions on Graphics, vol. 31, no. 4, pp. 107:1–107:9, 2012.
 [30] R. Weller and G. Zachmann, “Inner sphere trees for proximity and penetration queries.” in Robotics: Science and Systems, vol. 2, 2009.
 [31] B. Wang, F. Faure, and D. K. Pai, “Adaptive imagebased intersection volume,” in Proceedings of SIGGRAPH, 2012, pp. 97:1–97:9.
 [32] J.M. Lien, “Covering Minkowski sum boundary using points with applications,” Computer Aided Geometric Design, vol. 25, no. 8, pp. 652–666, 2008.
 [33] G. Nawratil, H. Pottmann, and B. Ravani, “Generalized penetration depth computation based on kinematical geometry,” Computer Aided Geometric Design, vol. 26, no. 4, pp. 425–443, May 2009.
 [34] L. Zhang, Y. J. Kim, and D. Manocha, “A fast and practical algorithm for generalized penetration depth computation,” in Robotics: Science and Systems, 2007.
 [35] Y. Kim, D. Manocha, and Y. Kim, “Hybrid penetration depth computation using local projection and machine learning,” in IEEE/RSJ International Conference on Intelligent Robots and Systems, 2015.
 [36] K. Kazerounian and J. Rastegar, “Object norms: A class of coordinate and metric independent norms for displacement,” in ASME Design Technical Conference, 1992, pp. 271–275.