A Large-scale Benchmark and an Inclusion-based Algorithm for Continuous Collision Detection A Large-scale Benchmark and an Inclusion-based Algorithm for Continuous Collision Detection

We crafted two datasets to compare the performance and correctness of CCD algorithms: (1) a handcrafted dataset that contains over 12 thousand point-triangle and 15 thousand edge-edge queries, and (2) a simulation dataset that contains over 18 million point-triangle and 41 million edge-edge queries. To foster replicability, we describe the format of the dataset in Appendix A.

The handcrafted queries are the union of queries simulated with Li et al. 2020 from the scenes in Erleben 2018 (Figure 2) and a set of handcrafted pairs for degenerate geometric configurations. These include: point-point degeneracies, near collisions (within a floating-point epsilon from collision), coplanar vertex-face and edge-edge motion (where the function $f$ (3) has infinite roots), degenerated function $F_{\text{vf}}$ and $F_{\text{ee}}$, and CCD queries with two or three roots.

The simulation queries were generated by running four nonlinear elasticity simulations. The first two simulations (Figure 3, top row) use the constraints of Verschoor and Jalba 2019 to simulate two cow heads colliding and a chain of rings falling. The second two simulations (Figure 3, bottom row) use the method of Li et al. 2020 to simulate a coarse mat twisting and the high-speed impact of a golf ball hitting a planar wall.

We compare seven state-of-the-art methods: (1) the interval root-finder (IRF) Snyder 1992, (2) the univariate interval root-finder (UIRF) (a special case of the rigid-body CCD from Redon et al. 2002), (3) the floating-point time-of-impact root finder (FPRF) Provot 1997 implemented in Vouga et al. 2010, (4) TightCCD (TCCD) Wang et al. 2015, (5) Root Parity (RP) Brochu et al. 2012, (6) a rational implementation of Root Parity (RRP) with the degenerate cases properly handled, and (7) BSC Tang et al. 2014. For each method, we collect the average query time, the number of false positives (i.e., there is no collision but the method detects one), and the number of false negatives (i.e., there is a collision but the method misses it). To obtain the ground truth, we solve the multivariate CCD formulation (Equations (1) and (2)) symbolically using Mathematica Wolfram Research Inc. 2020, which takes multiple seconds per query. Table 1 summarizes the results. Note that “Ours” corresponds to our new method that will be introduced and discussed in Section 5 and minimum separation floating-point time-of-impact root finder (MSRF) is a minimum separation CCD discussed in Section 6.2.

IRF. The inclusion-based root-finding described in Snyder 1992 can be applied to both the multivariate and univariate CCD. For the multivariate case, we can simply initialize the parameters of $F$ (i.e., $t, u, v$) with the size of the domain $\Omega$, evaluate $F$ and check if the origin is contained in the output interval Snyder et al. 1993. If it is, then we sequentially subdivide the parameters (thus shrinking the size of the intervals of $F$) until a user-tolerance $\delta$ is reached. In our comparison we use $\delta =10^{-6}$. The major advantage of this approach is that it is guaranteed to be conservative: It is impossible to shrink the interval of $F$ to zero. A second advantage is that a user can easily trade accuracy (number of false positives) for efficiency by simply increasing the tolerance $\delta$ (Appendix D). The main drawback is that bisecting $\Omega$ in the three dimensions makes the algorithm slow, and the use of interval arithmetic further increases the computational cost and prevents the use of certain compiler optimization techniques (such as instruction reordering). We implement this approach using the numerical type provided by the Boost interval library Schling 2011.

UIRF. Snyder 1992 can also be applied to the univariate function in Equation (3) by using the same subdivision technique on the single variable $t$ (as in Redon et al. 2002 but for linear trajectories). The result of this step is an interval containing the earliest root in $t$, which is then plugged inside a geometric predicate to check if the primitives intersect in that interval. While finding the roots with this approach might, at a first glance, seem easier than in the multi-variate case and thus more efficient, this is not the case in our experiments. If the polynomial has infinite roots, then this algorithm will have to refine the entire domain to the maximal allowed resolution, and check the validity of each interval, making it correct but very slow on degenerate cases (Appendix D). This results in a longer average runtime than its multivariate counterpart. Additionally, it is impossible to control the accuracy of the other two parameters (i.e., $u, v$), thus introducing more false positives.

FPRF. Vouga et al. 2010 aim to solve the univariate CCD problem using only floating-point computation. To mitigate false negatives, the method uses a numerical tolerance $\eta$ (Appendix E) shows how $\eta$ affects running time, the false positive, and negative). The major limitations are that the number of false positives cannot be directly controlled as it depends on the relative position of the input primitives and that false negatives can appear if the parameter is not tuned accordingly to the objects velocity and scale. Additionally, the reference implementation does not handle the edge-edge CCD when the two edges are parallel. This method is one of the fastest, which makes it a very popular choice in many simulation codes.

TCCD. TightCCD is a conservative floating-based implementation of Tang et al. 2014. It uses the univariate formulation coupled with three inequality constraints (two for the edge-edge case) to ensure that the univariate root is a CCD root. The algorithm expresses the cubic polynomial $f$ as a product and sum of three low-order polynomials in Bernstein form. With this reformulation the CCD problem becomes checking if univariate Bernstein polynomials are positive, which can be done by checking some specific points. This algorithm is extremely fast but introduces many false positives that are impossible to control. In our benchmark, this is the only non-interval method without false negatives. The major limitation of this algorithm is that it always detects collision if the primitives are moving in the same plane, independently from their relative position.

RP and RRP. These two methods use the multivariate formulation $F$ (Equations (1) and (2)). The main idea is that the parity of the roots of $F$ can be reduced to a ray casting problem. Let $\partial \Omega$ be the boundary of $\Omega$, the algorithm shoots a ray from the origin and counts the parity of the intersection between the ray and $F(\partial \Omega)$ that corresponds to the parity of the roots of $F$. Parity is however insufficient for CCD: These algorithms cannot differentiate between zero roots (no collision) and two roots (collision), since they have the same parity. We note that this is a rare case happening only with sufficiently large timesteps and/or velocities: We found 13 (handcrafted dataset) and 7 (simulation dataset) queries where these methods report a false negative.

We note that the algorithm described in Brochu et al. 2012 (and its reference implementation) does not handle some degenerate cases leading to both false negatives and positives. For instance, in Appendix B, we show an example of a “hourglass” configuration where RP misses the collision, generating a false negative. To overcome this limitations and provide a fair comparison to these techniques, we implemented a naïve version of this algorithm that handles all the degenerate cases using rational numbers to simplify the coding (see the additional materials). We opted for this rational implementation since properly handling the degeneracies using floating-point requires designing custom higher precision predicates for all cases. The main advantage of this method is that it is exact (when the degenerate cases are handled) as it does not contain any tolerance and thus has zero false positives. We note that the runtime of our rational implementation is extremely high and not representative of the runtime of a proper floating point implementation of this algorithm.

BSC. This efficient and exact method uses the univariate formulation coupled with inequality constraints to ensure that the coplanar primitives intersects. The coplanarity problem reduces to checking if $f$ in Bernstein form has a root. Tang et al. 2014 explain how this can be done exactly by classifying the signs of the four coefficients of the cubic Bernstein polynomial. The classification holds only if the cubic polynomial has monotone curvature; which can be achieved by splitting the curve at the inflection point. This splitting, however, cannot be computed exactly as it requires divisions (Appendix C). In our comparison, we modified the reference implementation to fix a minor typo in the code and to handle $f$ with inflection points by conservatively reporting collision. This change introduces potential false positives, and we refer to the additional material for more details and for the patch we applied to the code.

Discussion and Conclusions. From our extensive benchmark of CCD algorithms, we observe that most algorithms using the univariate formulation have false negatives. While the reduction to univariate root findings provides a performance boost, filtering the roots (without introducing false positives) is a challenging problem for which a robust solution is still elusive.

Surprisingly, only the oldest method, IRF, is at the same time reasonably efficient (e.g., it does not take multiple seconds per query as Mathematica), correct (i.e., no false negatives), and returns a small number of false positives (which can be controlled by changing the tolerance $\delta$). It is, however, slower than other state-of-the-art methods, which is likely the reason why it is currently not widely used. In the next section, we show that it is possible to change the inclusion function used by this algorithm to keep its favorable properties, while decreasing its runtime by ${\sim }250$ times, making its performance competitive with state-of-the-art methods.

An interval $i = [a, b]$ is defined as

An interval can be used to define an inclusion function. Formally, given an $m$-dimensional interval $D$ and a continuous function , an inclusion function for $g$, written , is a function such that

A convergent inclusion function can be used to find a root of a function $g$ over a domain bounded by the interval $I_0=[\mathcal {L}(x_1),\mathcal {R}(x_1)]\times \cdots \times [\mathcal {L}(x_m),\mathcal {R}(x_m)]$. To find the roots of $g$, we sequentially bisect the initial $m$-dimensional interval $I_0$, until it becomes sufficiently small (Algorithm 1). Figure 4 shows a one-dimensional (1D) example (i.e., ) of a bisection algorithm. The algorithm starts by initializing a stack $S$ of intervals to be checked with $I_0$ (line 3). At every level $\ell$ (line 5), the algorithm retrieves an interval $I$ from $S$ and evaluates the inclusion function to obtain the interval $I_g$ (line 7). Then it checks if the root is included in $I_g$ (line 8). If not, then $I$ can be safely discarded since $I_g$ bounds the range of $g$ over the domain bounded by $I$. Otherwise ($0 \in I_g$), it checks if $w(I)$ is smaller than a user-defined threshold $\delta$. If so, then it appends $I$ to the result (line 10). If $I$ is too large, then the algorithm splits one of its dimensions (e.g., $[\mathcal {L}(x_1),\mathcal {R}(x_1)]$ is split in $[\mathcal {L}(x_1),\tilde{x}_1]$ and $[\tilde{x}_1, \mathcal {R}(x_1)]$ with $\tilde{x}_1 = (\mathcal {L}(x_1)+\mathcal {R}(x_1))/2$) and appends the two new intervals $I_1, I_2$ to the stack $S$ (line 13).

Generic Construction of Inclusion Functions. Snyder 1992 proposes the use of interval arithmetic as a universal and automatic way to build inclusion functions for arbitrary expressions. However, interval arithmetic adds a performance overhead to the computation. For example, the product between two intervals is

Instead of using interval arithmetic to construct the inclusion function for the interval $ I_\Omega = I_t \times I_u \times I_v = [0, 1] \times [0, 1] \times [0, 1]$ around the domain $\Omega$, we propose to define an inclusion function tailored for $F$ (both for Equation (1) and (2)) as the box

The inclusion function $B_F$ turns out to be ideal for constructing a predicate: to use this inclusion function in the solve algorithm (Algorithm 1), we only need to check if, for a given interval $I$, $B_F(I)$ contains the origin (line 8). Such a Boolean predicate can be conservatively evaluated using floating point filtering.

Conservative Predicate Evaluation. Checking if the origin is contained in an axis-aligned box is trivial and it reduces to checking if the zero is contained in the three intervals defining the sides of the box. In our case, this requires us to evaluate the sign of $F$ at the eight box corners. However, the vertices of the co-domain are computed using floating point arithmetic and can thus be inaccurate. We use forward error analysis to conservatively account for these errors as follows.

Without loss of generality, we focus only on the $x$-axis. Let $\lbrace v_i^x\rbrace , {i=1,\ldots , 8}$ be the set of $x$-coordinates of the eight vertices of the box represented in double precision floating-point numbers. The error bound for $F$ (on the $x$-axis) is

Efficient Evaluation. The $x,y,z$ predicates defined above depend only on a subset of the coordinates of the eight corners of $B_F(I)$. We can optimally vectorize the evaluation of the eight corners using AVX2 instructions (${\sim }4\times$ improvement in performance), since it needs to be evaluated on eight points and all the computation is standard floating-point arithmetic. Note that we used AVX2 instructions because newer versions still have spotty support on current processors. After the eight points are evaluated in parallel, applying the floating-point filter involves only a few comparisons. To further reduce computation, we check one axis at a time and immediately return if any of the intervals do not contain the origin.

Algorithm.

We describe our complete algorithm in pseudocode in Algorithm 2. The input to our algorithm are the eight points representing two primitives (either vertex-face or edge-edge), a user-controlled numerical tolerance $\delta \gt 0$ (if not specified otherwise, in the experiment we use the default value $\delta =10^{-6}$), and the maximum number of checks $m_I \gt 0$ (we use the default value $m_I=10^6$). These choice are based on our empirical results (figures 8 and 9). The output is a conservative estimate of the earliest time of impact or infinity if the two primitives do not collide in the time intervals coupled with the reached tolerance.

Our algorithm iteratively checks the box $B=B_F(I)$, with $I=I_t\times I_u \times I_v$ = $[t_1, t_2]\times [u_1, u_2]\times [v_1, v_2]\subset I_\Omega$ (initialized with $[0,1]^3$). To guarantee a uniform box size while allowing early termination of the algorithm, we explore the space in a breadth-first manner and record the current explored level $\ell$ (line 6). Since our algorithm is designed to find the earliest time of impact, we sort the visiting queue $Q$ with respect to time (line 21).

At every iteration, we check if $B$ intersects the cube $C_{\varepsilon } = [-\varepsilon ^x, \varepsilon ^x]\times [-\varepsilon ^y, \varepsilon ^y]\times [-\varepsilon ^z, \varepsilon ^z]$ (line 9); if it does not, then we can safely ignore $I$ since there are no collisions.

If $B \cap C_{\varepsilon } \ne \emptyset$, then we first check if $w(B) \lt \delta$ or if $B$ is contained inside the $\varepsilon$-box (line 15). In this case, it is unnecessary to refine the interval $I$ more since it is either already small enough (if $w(B) \lt \delta$) or any refinement will lead to collisions (if $B \subseteq C_{\varepsilon }$). We return $I_t^l$ (i.e., the left hand-side of the $t$ interval of $I$) only if $I$ was the first intersecting interval of this current level (line 16). If $I$ is not the first intersecting in the current level, then there is an intersecting box (which is larger than $\delta$) with an earlier time since the queue is sorted according to time (Figure 5).

If $B$ is too big, then we split the interval $I$ in two sub-intervals and push them to the priority queue $Q$ (line 19). Note that, differently from Algorithm 1, we use a priority queue $Q$ instead of the stack $S$. For the vertex-triangle CCD, the domain $\Omega$ is a prism, thus, after spitting the interval (line 19), we append $I_1, I_2$ to $Q$ only if they intersect with $\Omega$. To ensure that $B$ shrinks uniformly (since the termination criteria, Line 15, is $w(B) \lt \delta$) we conservatively estimate the width of $B$ (in the codomain) from the widths of the domain's (i.e., where the algorithm is acting) intervals $I_t, I_u, I_v$:

Proof. While $B_F(I)$ is an interval, for the purpose of the proof we equivalently define it as an axis-aligned bounding box whose eight vertices are $b_i$. We will use the super-script notation to refer to the $x,y,z$ component of a 3D point (e.g., $b_i^x$ is the $x$-component of $b_i$) and define the set $\mathcal {I}=\lbrace 1,\ldots ,8\rbrace$. By using the box definition the width of $B_F(I)$ can be written as

\begin{equation*} w(B_F(I)) = \Vert b_M-b_m\Vert _\infty \end{equation*}

with

\begin{equation*} b_M^k=\max _{i\in \mathcal {I}}(b_i^k)\quad \text{and}\quad b_m^k=\min _{i\in \mathcal {I}}(b_i^k). \end{equation*}

Since $B_F(I)$ is the tightest axis-aligned inclusion function (Proposition 1)

\begin{equation*} b_M^k \leqslant \max _{i\in \mathcal {I}}{v_i^k}, \quad b_m^k \leqslant \min _{i\in \mathcal {I}}{v_i^k}, \end{equation*}

where $v_i = F(I_t^j, I_u^k, I_v^l),$ with $j,k,l\in \lbrace l,r\rbrace$, thus for any coordinate $k$ we bound

\begin{equation*} b_M^k-b_m^k= \max _{i,j\in \mathcal {I}}(v_i^k-v_j^k)\leqslant \max _{i,j\in \mathcal {I}}\Vert v_i-v_j\Vert _\infty . \end{equation*}

For any pair of $v_i$ and $v_j$ we have

\begin{equation*} v_i-v_j = s_1\alpha _{l,m}+s_2\beta _{n,p}+s_3\gamma _{p,q}, \end{equation*}

for some indices $l,m,n,o,p,q\in \lbrace 1,2\rbrace$ and constant $s_1,s_2,s_3 \in \lbrace -1,0,1\rbrace$ with

\begin{equation*} \alpha _{i,j}=w(I_t)\big (F(0,u_i,v_j)-F(1,u_i,v_j)\big), \end{equation*}

\begin{equation*} \beta _{i,j}=w(I_u)\big (F(t_i,0,v_j)-F(t_i,1,v_j)\big), \end{equation*}

\begin{equation*} \gamma _{i,j}=w(I_v)\big (F(t_i,u_j,0)-F(t_i,u_j,1)\big), \end{equation*}

since $F$ is linear on the edges. We note that $\alpha _{i,j}, \beta _{i,j}$, and $\gamma _{i,j}$ are the 12 edges of the box $B_F$. We now define

\begin{equation*} e_t^k = \max _{i,j\in \lbrace 1,2\rbrace } |\alpha ^k_{i,j}|,\quad e_u^k = \max _{i,j\in \lbrace 1,2\rbrace } |\beta ^k_{i,j}|,\quad e_v^k = \max _{i,j\in \lbrace 1,2\rbrace } |\gamma ^k_{i,j}| \end{equation*}

which allows us to bound

\begin{equation*} \max _{i,j\in \mathcal {I}}\Vert v_i-v_j\Vert _\infty \leqslant \Vert e_t+e_u+e_v\Vert _\infty \leqslant \Vert e_t\Vert _\infty +\Vert e_u\Vert _\infty +\Vert e_v\Vert _\infty . \end{equation*}

Since

\begin{equation*} \Vert e_t\Vert _\infty \leqslant w(I_t)\max _{i,j=1,2}\Vert F(t_1,u_i,v_j)-F(t_2,u_i,v_j)\Vert _\infty =w(I_t)\kappa _t/3, \end{equation*}

and similarly , we have

\begin{equation*} \Vert e_t\Vert _\infty +\Vert e_u\Vert _\infty +\Vert e_v\Vert _\infty \leqslant \frac{w(I_t)\kappa _t+w(I_u)\kappa _u+w(I_v)\kappa _v}{3} \end{equation*}

Finally, from the assumption ( 6) it follows that

\begin{equation*} w(B_F(I)) \leqslant \max _{i,j\in \mathcal {I}}\Vert v_i-v_j\Vert _\infty \leqslant \Vert e_t\Vert _\infty +\Vert e_u\Vert _\infty +\Vert e_v\Vert _\infty \lt \alpha . \end{equation*}

□

Using the estimate of the width of $I_t, I_u, I_v$ we split the dimension that leads to the largest estimated dimension in the range of $F$ (line 28).

Fixed Runtime or Fixed Accuracy. To ensure a bounded runtime, which is very useful in many simulation applications, we stop the algorithm after an user-controlled number of checks $m_I$. To ensure that our algorithm always returns a conservative time of impact we record the first colliding interval $I_f$ of every level (line 11). When the maximum number of check is reached we can safely return the latest recorded interval $I_f$ (line 13) (Figure 6). We note that our algorithm will not respect the user specified accuracy when it terminates early: If a constant accuracy is required by applications, then this additional termination criteria could be disabled, obtaining an algorithm with guaranteed accuracy but sacrificing the bound on the maximal running time. Note that without the termination criteria $m_I$, it is possible (while rare in our experiments) that the algorithm will take a long time to terminate, or run out of memory due to storing the potentially large list of candidate intervals $L$.

Our algorithm is implemented in C++ and uses Eigen Guennebaud et al. 2010 for the linear algebra routines (with the -avx2 g++ flag). We run our experiments on a 2.35 GHz AMD EPYC™ 7452. We attach the reference implementation and the data used for our experiments, which will be released publicly.

The running time of our method is comparable to the floating-point methods, while being provably correct, for any choice of parameters. For this comparison we use a default tolerance $\delta =10^{-6}$ and default number of iterations $m_I=10^6$. All queries in the simulation dataset terminate within $10^6$ checks, while for the handcrafted dataset only $0.25\%$ and $0.55\%$ of the vertex-face and edge-edge queries required more than $10^6$ checks, reaching an actual maximal tolerance $\delta$ of $2.14\times 10^{-5}$ and $6.41\times 10^{-5}$ for vertex-face and edge-edge respectively. We note that, despite the percentages begin small, by removing $m_I$ the handcrafted queries take 0.015774 and 0.042477 s on average for vertex-face and edge-edge respectively. This is due to the large number of degenerate queries, as can be seen from the long tail in the histogram of the run-times (Figure 7). We did not observe any noticeable change of running time for the simulation dataset.

Our algorithm has two user-controlled parameters ($\delta$ and $m_I$) to control the accuracy and running time. The tolerance $\delta$ provides a direct control on the achieved accuracy and provides an indirect effect on the running time (Figure 8). The other parameter, $m_I$, directly controls the maximal running time of each query: for small $m_I$ our algorithm will terminate earlier, resulting in a lower accuracy and thus more chances of false positives (Figure 9, top). We remark that, in practice, very few queries require so many subdivisions: by reducing $m_I$ to the very low value of 100, our algorithm early-terminates only on ${\sim }0.07$% of the 60 million queries in the simulation dataset.

The input to our MSCCD algorithm are the same as the standard CCD (eight coordinates, $\delta$, and $m_I$) and the minimum separation distance . Our algorithm returns the earliest time of impact indicating if two primitives become closer than $d$ as measured by the $L^\infty$ norm.

We wish to check whether the box $B_F(\Omega)$ intersects a cube of side $2d$ centered on the origin (Figure 10). Equivalently, we can construct another box $B_F^{\prime }(\Omega)$ by displacing the six faces of $B_F(\Omega)$ outward at a distance $d$, and then check whether this enlarged box contains the origin. This check can be done as for the standard CCD (Section 5), but the floating point filters must be recalculated to account for the additional sum (indeed, we add/subtract $d$ to/from all the coordinates). Hence, the filters for $F^{\prime }$ are as follows:

To account for minimum separations, the only change in our algorithm is at line 7 where we need to enlarge $B$ by $d$ and in lines 9 and 15 since $C_{\varepsilon }$ needs to be replaced with $C_\epsilon = [-\epsilon ^x, \epsilon ^x]\times [-\epsilon ^y, \epsilon ^y]\times [-\epsilon ^z, \epsilon ^z]$.

To the best of our knowledge, the MSRF Harmon et al. 2011 implemented in Lu et al. 2019, is the only public code supporting minimal separation queries. While not explicitly constructed for MSCCD, FPRF uses a distance tolerance to limit false negatives, similarly to an explicit minimum separation. We compare the results and performance in Appendix E. MSRF. Uses the univariate formulation, which requires to find the roots of a high-order polynomial, and it is thus unstable when implemented using floating-point arithmetic.

Table 2 reports timings, false positive, and false negatives for different separation distances $d$. As $d$ shrinks (around $10^{-16}$) the results of our method with MSCDD coincide with the ones with $d=0$ since the separation is small. For these small tolerances, MSRF runs into numerical problems and the number of false negatives increases. Figure 11 shows the average query time versus the separation distance $d$ for the simulation dataset, since our method only requires to check the intersection between boxes, the running time largely depends on the number of detected collision, and the average is only mildly affected by the choice of $d$.

In the traditional constraint based collision handling (such as that of Verschoor and Jalba 2019), collision response is handled by performing an implicit timestep as a constrained optimization. The goal is to minimize a elastic potential while avoiding interpenetration through gap constraints. To avoid handling all possible collisions during a simulation, a subset of active collisions constraints $C_A$ is usually constructed. This set not only avoids infeasibilities, but also improves performance by having fewer constraints. There are many activation strategies, but for the sake of brevity we focus here on the strategies used by Verschoor and Jalba 2019.

Algorithm 3 shows how CCD is used to compute the active set $C_A$. Given the starting and ending vertex positions, $x_0$ and $x_1$, we compute the time of impact for each collision candidate $c \in C$. We use the notation $x_i \cap c$ to indicate selecting the constrained vertices from $x_i$. If the candidate $c$ is an actual collision, that is, , then we add this constraint and the time of impact, $t$, to the active set, $C_A$.

From the active constraint set the constraints of Verschoor and Jalba 2019 are computed as

Because of the difficulty for a simulation solver to maintain and not violate constraints, it is common to offset the constraints such that

Figure 12 shows example of simulations run with different numerical tolerance $\delta$. Changing $\delta$ has little effect on the simulation in terms of run-time, but for large values of $\delta$, it can affect accuracy. We observe that for a the simulation is more likely to contain intersections. This is most likely due to the inaccuracies in the contact points used in the constraints.

A line search is used in a optimization to ensure that every update decreases the energy $E$. That is, given an update, $\Delta x$, to the optimization variable $x$, we want to find a step size $\alpha$ such that $E(x + \alpha \Delta x) \lt E(x)$. This ensure that we make progress toward a minimum.

When used in a line search algorithm, CCD can be used to prevent intersections and tunneling. This requires modifying the maximum step length to the time of impact. As observed by Li et al. 2020, the standard CCD formulation without minimal separation cannot be used directly in a line search algorithm. Let $t^\star$ the earliest time of impact (i.e., $F(t^\star , \tilde{u}, \tilde{v}) = 0$ for some $\tilde{u}, \tilde{v}$ and there is no collision between 0 and $t^\star$) and assume that the energy at $E(x_0 + t^\star \Delta x) \lt E(x_0)$ (Algorithm 4, line 22). In this case the step $\alpha = t^\star$ is a valid descent step that will be used to update the position $x$ in outer iteration (e.g., Newton optimization loop). In the next iteration, the line search will be called with the updated position and the earliest time of impact will be zero since we selected $t^\star$ in the previous iteration. This prevents the optimization from making progress because any direction $\Delta x$ will lead to a time of impact $t=0$. To avoid this problem we need the line search to find an appropriate step-size $\alpha$ along the update direction that leaves “sufficient space” for the next iteration, so that the barrier in Li et al. 2020 will be active and steer the optimization away from the contact position. Formally, we aim at finding a valid CCD sequence $\lbrace t_i\rbrace$ such that

Constructing a Sequence. Let $0\lt p\lt 1$ be a user-defined tolerance ($p$ close to 1 will produce a sequence $\lbrace t_i\rbrace$ converging faster) and $d_i$ be the distance between two primitives. We propose to set $d=p d_i$, and ensure that no primitive are closer than $d$. Without loss of generality, we assume that $F(x+\Delta x) =0$, that is, taking the full step will lead to contact. By taking successive steps in the same direction, $d_i$ will shrink to zero ensuring $t_i$ to converge to $t^\star$. Similarly we will obtain a growing sequence $t_i$ since $d$ decreases as we proceed with the iterations. Finally, it is easy to see that $p=t_i/t_{i+1}$, which can be close to 1.

To account for the aforementioned problem, we propose to use our MSCCD algorithm to return a valid CCD sequence when employed in a line search scenario. For a step $i$, we define $\delta ^i$ as the tolerance, $\epsilon _i$ the numerical error (8), and $\rho _i$ as the maximum numerical error in computing the distances $d_i$ from the candidates set $C$ (line 5). $\rho _i$ should be computed using forward error analysis on the application-specific distance computation: Since the applications are not the focus of our article, we used a fixed $\rho _i=10^{-9}$, and we leave the complete forward analysis as a future work. (We note that our approximation might thus introduce zero length steps, this however did not happen in our experiments.) If $d_i -(\delta _i +\epsilon _i + \rho _i) \gt d$, then our MSCCD is guaranteed to find a time of impact larger than zero. Thus, if we set $d=p d_i$ (line 7), then we are guaranteed to find a positive time of impact if

As visible from Table 2, our MSCCD slows down as $d$ grows. Since the actual minimum distance is not relevant in the line search algorithm, our experiments suggest to cap it at $\delta$ (line 7). To avoid unnecessary computations and speedup the MSCCD computations, our algorithm, as suggested by Redon et al. 2002, can be easily modified to accept a shorter time interval (line 13): It only requires to change the initialization of $I$ (Algorithm 2 line 3). These two modifications lead to a $8\times$ speedup in our experiments. We refer to this algorithm with MSCCD (i.e., Algorithm 2 with MSCDD, Section 6.1, and modified initialization of $I$) as SolveMSCCD.

Figure 13 shows a simulation using our MSCCD in line search to keep the bodies from intersecting for different $\delta$. As illustrated in the previous section, the effect of $\delta$ is negligible as long as . Timings vary depending on the maximum number of iterations. Because the distance $d$ varies throughout the simulation, some steps take longer than others (as seen in Figure 11). We note that if we use the standard CCD formulation $F=0$, then the line search gets stuck in all our experiments, and we were not able to find a solution. Note that for a line search based method it is crucial to have a conservative CCD/MSCCD algorithm: The videos in the additional material shows that a false negative leads to an artefact in the simulation.