forked from BachiLi/diffvg
-
Notifications
You must be signed in to change notification settings - Fork 0
/
solve.h
59 lines (55 loc) · 1.51 KB
/
solve.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
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;
}
}