Add max_iterations to additive CCD (#145)

Author: zfergus

python/src/ccd/additive_ccd.cpp

@@ -9,13 +9,14 @@ void define_additive_ccd(py::module_& m)
 {
     py::class_<AdditiveCCD, NarrowPhaseCCD>(m, "AdditiveCCD")
         .def(
-            py::init<const double>(),
+            py::init<const double, const double>(),
             R"ipc_Qu8mg5v7(
             Construct a new AdditiveCCD object.
 
             Parameters:
                 conservative_rescaling: The conservative rescaling of the time of impact.
             )ipc_Qu8mg5v7",
+            py::arg("max_iterations") = AdditiveCCD::DEFAULT_MAX_ITERATIONS,
             py::arg("conservative_rescaling") =
                 AdditiveCCD::DEFAULT_CONSERVATIVE_RESCALING)
         .def_readonly_static(

src/ipc/candidates/candidates.cpp

@@ -317,8 +317,8 @@ double Candidates::compute_cfl_stepsize(
     // If alpha_F < 0.5 * alpha_C, then we should do full CCD.
     if (alpha_F < 0.5 * alpha_C) {
         return ipc::compute_collision_free_stepsize(
-            mesh, vertices_t0, vertices_t1, min_distance, //
-            broad_phase_method, narrow_phase_ccd);
+            mesh, vertices_t0, vertices_t1, min_distance, broad_phase_method,
+            narrow_phase_ccd);
     }
     return std::min(alpha_C, alpha_F);
 }

src/ipc/ccd/additive_ccd.cpp

@@ -9,6 +9,7 @@
 //  • return true if the initial distance is less than the minimum distance
 //  • add an explicit tmax parameter rather than relying on the initial value of
 //    toi
+//  • add a maximum number of iterations to limit the computation time
 //
 // NOTE: These methods are provided for reference comparison with [Li et al.
 // 2021] and is not utilized by the high-level functionality. In compairson to
@@ -60,8 +61,10 @@ namespace {
     }
 } // namespace
 
-AdditiveCCD::AdditiveCCD(const double _conservative_rescaling)
-    : conservative_rescaling(_conservative_rescaling)
+AdditiveCCD::AdditiveCCD(
+    const long _max_iterations, const double _conservative_rescaling)
+    : max_iterations(_max_iterations)
+    , conservative_rescaling(_conservative_rescaling)
 {
 }
 
@@ -72,8 +75,7 @@ bool AdditiveCCD::additive_ccd(
     const double max_disp_mag,
     double& toi,
     const double min_distance,
-    const double tmax,
-    const double conservative_rescaling)
+    const double tmax) const
 {
     assert(conservative_rescaling > 0 && conservative_rescaling <= 1);
 
@@ -87,9 +89,14 @@ bool AdditiveCCD::additive_ccd(
     assert(d_func > 0);
     const double gap = // (d - ξ) = (d² - ξ²) / (d + ξ)
         (1 - conservative_rescaling) * d_func / (d + min_distance);
+    if (gap < std::numeric_limits<double>::epsilon()) {
+        logger().warn(
+            "Small gap {:g} ≤ ε in Additive CCD can lead to missed collisions",
+            gap);
+    }
 
     toi = 0;
-    while (true) {
+    for (long i = 0; max_iterations < 0 || i < max_iterations; ++i) {
         // tₗ = η ⋅ (d - ξ) / lₚ = η ⋅ (d² - ξ²) / (lₚ ⋅ (d + ξ))
         const double toi_lower_bound = conservative_rescaling * d_func
             / ((d + min_distance) * max_disp_mag);
@@ -108,6 +115,12 @@ bool AdditiveCCD::additive_ccd(
         if (toi > tmax) {
             return false; // collision occurs after tmax
         }
+
+        if (max_iterations < 0 && i == DEFAULT_MAX_ITERATIONS) {
+            logger().warn(
+                "Slow convergence in Additive CCD. Perhaps the gap is too small (gap={:g})?",
+                gap);
+        }
     }
 
     return true;
@@ -151,8 +164,7 @@ bool AdditiveCCD::point_point_ccd(
     const VectorMax12d dx = stack(dp0, dp1);
 
     return additive_ccd(
-        x, dx, distance_squared, max_disp_mag, toi, min_distance, tmax,
-        conservative_rescaling);
+        x, dx, distance_squared, max_disp_mag, toi, min_distance, tmax);
 }
 
 bool AdditiveCCD::point_edge_ccd(
@@ -199,8 +211,7 @@ bool AdditiveCCD::point_edge_ccd(
     const VectorMax12d dx = stack(dp, de0, de1);
 
     return additive_ccd(
-        x, dx, distance_squared, max_disp_mag, toi, min_distance, tmax,
-        conservative_rescaling);
+        x, dx, distance_squared, max_disp_mag, toi, min_distance, tmax);
 }
 
 bool AdditiveCCD::point_triangle_ccd(
@@ -248,8 +259,7 @@ bool AdditiveCCD::point_triangle_ccd(
     const VectorMax12d dx = stack(dp, dt0, dt1, dt2);
 
     return additive_ccd(
-        x, dx, distance_squared, max_disp_mag, toi, min_distance, tmax,
-        conservative_rescaling);
+        x, dx, distance_squared, max_disp_mag, toi, min_distance, tmax);
 }
 
 bool AdditiveCCD::edge_edge_ccd(
@@ -310,8 +320,7 @@ bool AdditiveCCD::edge_edge_ccd(
     const VectorMax12d dx = stack(dea0, dea1, deb0, deb1);
 
     return additive_ccd(
-        x, dx, distance_squared, max_disp_mag, toi, min_distance, tmax,
-        conservative_rescaling);
+        x, dx, distance_squared, max_disp_mag, toi, min_distance, tmax);
 }
 
 } // namespace ipc

src/ipc/ccd/additive_ccd.hpp

@@ -17,13 +17,18 @@ namespace ipc {
 /// @brief Additive Continuous Collision Detection (CCD) from [Li et al. 2021].
 class AdditiveCCD : public NarrowPhaseCCD {
 public:
+    /// The default maximum number of iterations used with Tight-Inclusion CCD.
+    static constexpr long DEFAULT_MAX_ITERATIONS = 10'000'000l;
+    /// Unlimitted number of iterations.
+    static constexpr long UNLIMITTED_ITERATIONS = -1;
     /// The default conservative rescaling value used to avoid taking steps
     /// exactly to impact. Value choosen to based on [Li et al. 2021].
     static constexpr double DEFAULT_CONSERVATIVE_RESCALING = 0.9;
 
     /// @brief Construct a new AdditiveCCD object.
     /// @param conservative_rescaling The conservative rescaling of the time of impact.
     AdditiveCCD(
+        const long max_iterations = UNLIMITTED_ITERATIONS,
         const double conservative_rescaling = DEFAULT_CONSERVATIVE_RESCALING);
 
     /// @brief Computes the time of impact between two points using continuous collision detection.
@@ -122,6 +127,9 @@ class AdditiveCCD : public NarrowPhaseCCD {
         const double min_distance = 0.0,
         const double tmax = 1.0) const override;
 
+    /// @brief Maximum number of iterations.
+    long max_iterations;
+
     /// @brief The conservative rescaling value used to avoid taking steps exactly to impact.
     double conservative_rescaling;
 
@@ -133,15 +141,14 @@ class AdditiveCCD : public NarrowPhaseCCD {
     /// @param tmax The maximum time to check for collisions.
     /// @param conservative_rescaling The amount to rescale the objects by to ensure conservative advancement.
     /// @return True if a collision was detected, false otherwise.
-    static bool additive_ccd(
+    bool additive_ccd(
         VectorMax12d x,
         const VectorMax12d& dx,
         const std::function<double(const VectorMax12d&)>& distance_squared,
         const double max_disp_mag,
         double& toi,
         const double min_distance = 0.0,
-        const double tmax = 1.0,
-        const double conservative_rescaling = DEFAULT_CONSERVATIVE_RESCALING);
+        const double tmax = 1.0) const;
 };
 
 } // namespace ipc

tests/src/tests/ccd/test_point_edge_ccd.cpp

@@ -29,7 +29,8 @@ void check_toi(
     CHECK(is_colliding);
     CHECK(toi <= toi_expected);
 
-    const AdditiveCCD additive_ccd(/*conservative_rescaling=*/0.99);
+    const AdditiveCCD additive_ccd(
+        AdditiveCCD::UNLIMITTED_ITERATIONS, /*conservative_rescaling=*/0.99);
     is_colliding = additive_ccd.point_edge_ccd(
         p_t0, e0_t0, e1_t0, p_t1, e0_t1, e1_t1, toi);
     CHECK(is_colliding);
@@ -170,7 +171,8 @@ TEST_CASE("Point-edge CCD", "[ccd][point-edge]")
     //     CHECK(toi <= toi_expected);
     // }
 
-    const AdditiveCCD additive_ccd(/*conservative_rescaling=*/0.99);
+    const AdditiveCCD additive_ccd(
+        AdditiveCCD::DEFAULT_MAX_ITERATIONS, /*conservative_rescaling=*/0.99);
     is_colliding = additive_ccd.point_edge_ccd(
         p_t0, e0_t0, e1_t0, p_t1, e0_t1, e1_t1, toi);
     REQUIRE(is_colliding == is_collision_expected);

tests/src/tests/ccd/test_point_point_ccd.cpp

@@ -28,7 +28,8 @@ TEST_CASE("Point-point CCD", "[ccd][point-point]")
     CHECK(is_colliding);
     CHECK(toi == Catch::Approx(0.5).margin(1e-3));
 
-    const AdditiveCCD additive_ccd(/*conservative_rescaling=*/0.999);
+    const AdditiveCCD additive_ccd(
+        AdditiveCCD::UNLIMITTED_ITERATIONS, /*conservative_rescaling=*/0.999);
     is_colliding = additive_ccd.point_point_ccd(
         p0_t0, p1_t0, p0_t1, p1_t1, toi, min_distance);