massive refactor and add TOG impl reference

Author: ryichando

README.md

@@ -52,6 +52,7 @@ Published in [ACM Transactions on Graphics (TOG)](https://dl.acm.org/doi/abs/10.
 
 ## 📝 Change History
 
+- (2025.02.28) Added a [reference branch and a Docker image of our TOG paper](#-technical-materials).
 - (2025.2.26) Added Floating Point-Rounding Errors in ACCD in [hindsight](./articles/hindsight.md).
 - (2025.2.7) Updated the [trapped example](./examples/trapped.ipynb) [[Video]](https://drive.google.com/file/d/1Qek0e0qBNWPlBb1hSOZ6o_e2Cqf5rGst/view) with squishy balls.
 - (2025.1.8) Added a [domino example](./examples/domino.ipynb) [[Video]](https://drive.google.com/file/d/1N9y8eZrjSQhAUhKwiO9w8jW_T18zPnYf/view).
@@ -64,7 +65,6 @@ Published in [ACM Transactions on Graphics (TOG)](https://dl.acm.org/doi/abs/10.
 - (2024.12.18) Added a [hindsight](./articles/hindsight.md) noting that the tilt angle was not $30^\circ$, but rather $26.57^\circ$
 - (2024.12.16) Removed thrust dependencies to fix runtime errors for the driver version `560.94` [[Issue Link]](https://github.com/st-tech/ppf-contact-solver/issues/1)
 
-
 ## 🎓 Technical Materials
 
 - 🎥 Main video [[Video]](https://drive.google.com/file/d/1OzPbUoqddUYDvXMvRnUHH7kz0nZhmt7K/view?usp=drive_link)
@@ -75,6 +75,20 @@ Published in [ACM Transactions on Graphics (TOG)](https://dl.acm.org/doi/abs/10.
 - 🤖 Supplementary scripts [[Directory]](https://drive.google.com/drive/folders/13CO068xLkd6ZSxsqtJQdNadgMrbbfSug?usp=drive_link)
 - 🔍 Singular-value eigenanalysis [[Markdown]](./articles/eigensys.md)
 
+#### 📌 Reference Implementation
+
+The main branch is undergoing frequent updates and may introduce breaking 🚧 changes relative to the paper.
+✨ Although the latest version offers a superior experience, we have created a new branch, ```sigasia-2024```, to retain consistency with the paper.
+
+- 🛠️ Only maintenance updates are planned for this branch.
+- 🚫 All algorithmic changes listed in this [[Markdown]](./articles/bug.md) are excluded from this branch.
+
+- 📦 We also provide a pre-compiled Docker image: ```ghcr.io/st-tech/ppf-contact-solver-compiled-sigasia-2024:latest``` of this branch.
+
+- 🌐 [Template Link for vast.ai](https://cloud.vast.ai/?ref_id=85288&creator_id=85288&name=ppf-contact-solver-sigasia-2024).
+
+- 🌐 [Template Link for RunPods](https://runpod.io/console/deploy?template=ooqpniuixi&ref=bhy3csxy)
+
 ## ⚡️ Requirements
 
 - 🔥 A modern NVIDIA GPU (Turing or newer)

src/args.rs

@@ -165,6 +165,15 @@ pub struct Args {
     #[clap(long, default_value_t = 0.01)]
     pub ccd_reduction: f32,
 
+    // Name: Maximum Search Direction Velocity
+    // Recommended Range: 5 m/s to 50 m/s
+    // Description:
+    // This parameter defines the maximum allowable search direction magnitude
+    // during optimization processes. It helps in controlling the step size
+    // and ensuring stability.
+    #[clap(long, default_value_t = 10.0)]
+    pub max_search_dir_vel: f32,
+
     // Name: Epsilon for Eigenvalue Analysis
     // Recommended Range: 1e-3 to 1e-2
     // Description:

src/builder.rs

@@ -346,6 +346,7 @@ pub fn make_param(args: &Args) -> data::ParamSet {
         line_search_max_t: args.line_search_max_t,
         ccd_tol: args.contact_ghat,
         ccd_reduction: args.ccd_reduction,
+        max_search_dir_vel: args.max_search_dir_vel,
         eiganalysis_eps: args.eiganalysis_eps,
         friction: args.friction,
         friction_eps: args.friction_eps,

src/bvh.rs

@@ -89,14 +89,34 @@ fn subdivide(mut obj: Vec<(Aabb, usize)>, node: &mut Vec<Node>) -> Node {
         let mean = compute_mean(obj.iter().map(|(aabb, _)| aabb));
         let var = compute_var(obj.iter().map(|(aabb, _)| aabb), mean);
         let (dim, _) = var.argmax();
-        obj.sort_by(|a, b| {
+        let mut left_obj = Vec::new();
+        let mut right_obj = Vec::new();
+        let mut deferred = Vec::new();
+        obj.iter().for_each(|&e| {
+            if e.0.max[dim] < mean[dim] {
+                left_obj.push(e);
+            } else if e.0.min[dim] > mean[dim] {
+                right_obj.push(e);
+            } else {
+                deferred.push(e);
+            }
+        });
+        deferred.sort_by(|a, b| {
             let a = a.0.min[dim] + a.0.max[dim];
             let b = b.0.min[dim] + b.0.max[dim];
             a.partial_cmp(&b).unwrap()
         });
-        let obj_len = obj.len();
-        let left_obj = Vec::from_iter(obj.drain(..obj_len / 2));
-        let right_obj = obj;
+        let deffered_len = deferred.len();
+        if deffered_len == 1 {
+            if left_obj.len() < right_obj.len() {
+                left_obj.push(deferred.pop().unwrap());
+            } else {
+                right_obj.push(deferred.pop().unwrap());
+            }
+        } else if deffered_len > 0 {
+            left_obj.extend(deferred.drain(..deffered_len / 2));
+            right_obj.extend(deferred);
+        }
         let (left_node, right_node) = (subdivide(left_obj, node), subdivide(right_obj, node));
         let left_index = node.len();
         node.push(left_node);

src/cpp/Makefile

@@ -14,7 +14,7 @@ CPP_FLAGS=-O3 -std=c++1z -Wall -Wno-unused-function -Wno-unused-function -Wno-st
 EIGEN_PATH=/usr/include/eigen3
 EIGSYS_PATH=../../eigsys
 
-BASE_DIRS=main vec utility csrmat contact aabb energy eigenanalysis barrier strainlimiting solver
+BASE_DIRS=main vec utility csrmat contact energy eigenanalysis barrier strainlimiting solver
 CUDA_FILES=$(foreach dir,$(BASE_DIRS),$(dir)/$(dir).cu)
 CPP_FILES=$(foreach dir,$(BASE_DIRS),$(dir)/$(dir).cpp)
 SIMPLELOG_FILES=simplelog/SimpleLog.cpp

src/cpp/aabb/aabb.cu

@@ -0,0 +1,101 @@
+// File: aabb.hpp
+// Author: Ryoichi Ando ([email protected])
+// License: Apache v2.0
+
+#ifndef AABB_DEF_HPP
+#define AABB_DEF_HPP
+
+#include "../data.hpp"
+
+#define AABB_MAX_QUERY 128
+
+namespace aabb {
+
+__device__ AABB join(const AABB &a, const AABB &b) {
+    return {Vec3f(std::min(a.min[0], b.min[0]), std::min(a.min[1], b.min[1]),
+                  std::min(a.min[2], b.min[2])),
+            Vec3f(std::max(a.max[0], b.max[0]), std::max(a.max[1], b.max[1]),
+                  std::max(a.max[2], b.max[2]))};
+}
+
+__device__ AABB make(const Vec3f &x0, const Vec3f &x1, const Vec3f &x2,
+                     float margin) {
+    AABB result = {Vec3f(std::min(x0[0], std::min(x1[0], x2[0])),
+                         std::min(x0[1], std::min(x1[1], x2[1])),
+                         std::min(x0[2], std::min(x1[2], x2[2]))),
+                   Vec3f(std::max(x0[0], std::max(x1[0], x2[0])),
+                         std::max(x0[1], std::max(x1[1], x2[1])),
+                         std::max(x0[2], std::max(x1[2], x2[2])))};
+    if (margin) {
+        for (unsigned i = 0; i < 3; ++i) {
+            result.min[i] -= margin;
+            result.max[i] += margin;
+        }
+    }
+    return result;
+}
+
+__device__ AABB make(const Vec3f &x0, const Vec3f &x1, float margin) {
+    AABB result = {Vec3f(std::min(x0[0], x1[0]), std::min(x0[1], x1[1]),
+                         std::min(x0[2], x1[2])),
+                   Vec3f(std::max(x0[0], x1[0]), std::max(x0[1], x1[1]),
+                         std::max(x0[2], x1[2]))};
+    if (margin) {
+        for (unsigned i = 0; i < 3; ++i) {
+            result.min[i] -= margin;
+            result.max[i] += margin;
+        }
+    }
+    return result;
+}
+
+__device__ AABB make(const Vec3f &x, float margin) {
+    return {
+        x - Vec3f(margin, margin, margin),
+        x + Vec3f(margin, margin, margin),
+    };
+}
+
+__device__ bool overlap(const AABB &a, const AABB &b) {
+    return (a.min[0] <= b.max[0] && a.max[0] >= b.min[0]) &&
+           (a.min[1] <= b.max[1] && a.max[1] >= b.min[1]) &&
+           (a.min[2] <= b.max[2] && a.max[2] >= b.min[2]);
+}
+
+template <typename F, typename T>
+__device__ static unsigned query(const BVH &bvh, Vec<AABB> aabb, F op,
+                                 T query) {
+    unsigned stack[AABB_MAX_QUERY];
+    unsigned count = 0;
+    unsigned head = 0;
+    if (bvh.node.size) {
+        stack[head++] = bvh.node.size - 1;
+        while (head) {
+            unsigned index = stack[--head];
+            if (op.test(aabb[index], query)) {
+                if (bvh.node[index][1] == 0) {
+                    unsigned leaf_index = bvh.node[index][0] - 1;
+                    if (op.test(aabb[index], query)) {
+                        if (op(leaf_index)) {
+                            count++;
+                        }
+                    }
+                } else {
+                    if (head + 2 >= AABB_MAX_QUERY) {
+                        printf("stack overflow!\n");
+                        assert(false);
+                        break;
+                    } else {
+                        stack[head++] = bvh.node[index][0] - 1;
+                        stack[head++] = bvh.node[index][1] - 1;
+                    }
+                }
+            }
+        }
+    }
+    return count;
+}
+
+} // namespace aabb
+
+#endif

src/cpp/aabb/aabb.hpp

@@ -5,13 +5,13 @@
 #ifndef ACCD_HPP
 #define ACCD_HPP
 
-#include "../math/distance.hpp"
+#include "distance.hpp"
 #include <cassert>
 
 namespace accd {
 
-template <class T> __device__ void centerize(SMat<T, 3, 4> &x) {
-    SVec<T, 3> mov = SVec<T, 3>::Zero();
+template <class T> __device__ void centerize(Mat3x4<T> &x) {
+    Vec3f mov = Vec3f::Zero();
     T scale(0.25f);
     for (int k = 0; k < 4; k++) {
         mov += scale * x.col(k);