60 lines
1.5 KiB
C++
60 lines
1.5 KiB
C++
#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;
|
|
}
|
|
}
|