initial commit

2020-09-03 22:30:30 -04:00
commit 413a3e5cee
148 changed files with 138536 additions and 0 deletions
--- a/solve.h
+++ b/solve.h
@@ -0,0 +1,59 @@
+#pragma once
+
+#include "diffvg.h"
+
+template <typename T>
+DEVICE
+inline bool solve_quadratic(T a, T b, T c, T *t0, T *t1) {
+    // From https://github.com/mmp/pbrt-v3/blob/master/src/core/pbrt.h#L419
+    T discrim = square(b) - 4 * a * c;
+    if (discrim < 0) {
+        return false;
+    }
+    T root_discrim = sqrt(discrim);
+
+    T q;
+    if (b < 0) {
+        q = -0.5f * (b - root_discrim);
+    } else {
+        q = -0.5f * (b + root_discrim);
+    }
+    *t0 = q / a;
+    *t1 = c / q;
+    if (*t0 > *t1) {
+        swap_(*t0, *t1);
+    }
+    return true;
+}
+
+template <typename T>
+DEVICE
+inline int solve_cubic(T a, T b, T c, T d, T t[3]) {
+    if (fabs(a) < 1e-6f) {
+        if (solve_quadratic(b, c, d, &t[0], &t[1])) {
+            return 2;
+        } else {
+            return 0;
+        }
+    }
+    // normalize cubic equation
+    b /= a;
+    c /= a;
+    d /= a;
+    T Q = (b * b - 3 * c) / 9.f;
+    T R = (2 * b * b * b - 9 * b * c + 27 * d) / 54.f;
+    if (R * R < Q * Q * Q) {
+        // 3 real roots
+        T theta = acos(R / sqrt(Q * Q * Q));
+        t[0] = -2.f * sqrt(Q) * cos(theta / 3.f) - b / 3.f;
+        t[1] = -2.f * sqrt(Q) * cos((theta + 2.f * T(M_PI)) / 3.f) - b / 3.f;
+        t[2] = -2.f * sqrt(Q) * cos((theta - 2.f * T(M_PI)) / 3.f) - b / 3.f;
+        return 3;
+    } else {
+        T A = R > 0 ? -pow(R + sqrt(R * R - Q * Q * Q), T(1./3.)):
+                           pow(-R + sqrt(R * R - Q * Q * Q), T(1./3.));
+        T B = fabs(A) > 1e-6f ? Q / A : T(0);
+        t[0] = (A + B) - b / T(3);
+        return 1;
+    }
+}