1 #include "geometry.hpp"
11 const glm::dmat4 &b_m,
15 glm::dvec3 a_corners[8] = {
16 glm::dvec3(a_m * glm::dvec4(a_box.min.x, a_box.min.y, a_box.min.z, 1.0)),
17 glm::dvec3(a_m * glm::dvec4(a_box.min.x, a_box.min.y, a_box.max.z, 1.0)),
18 glm::dvec3(a_m * glm::dvec4(a_box.min.x, a_box.max.y, a_box.min.z, 1.0)),
19 glm::dvec3(a_m * glm::dvec4(a_box.min.x, a_box.max.y, a_box.max.z, 1.0)),
20 glm::dvec3(a_m * glm::dvec4(a_box.max.x, a_box.min.y, a_box.min.z, 1.0)),
21 glm::dvec3(a_m * glm::dvec4(a_box.max.x, a_box.min.y, a_box.max.z, 1.0)),
22 glm::dvec3(a_m * glm::dvec4(a_box.max.x, a_box.max.y, a_box.min.z, 1.0)),
23 glm::dvec3(a_m * glm::dvec4(a_box.max.x, a_box.max.y, a_box.max.z, 1.0)),
26 glm::dvec3 b_corners[8] = {
27 glm::dvec3(b_m * glm::dvec4(b_box.min.x, b_box.min.y, b_box.min.z, 1.0)),
28 glm::dvec3(b_m * glm::dvec4(b_box.min.x, b_box.min.y, b_box.max.z, 1.0)),
29 glm::dvec3(b_m * glm::dvec4(b_box.min.x, b_box.max.y, b_box.min.z, 1.0)),
30 glm::dvec3(b_m * glm::dvec4(b_box.min.x, b_box.max.y, b_box.max.z, 1.0)),
31 glm::dvec3(b_m * glm::dvec4(b_box.max.x, b_box.min.y, b_box.min.z, 1.0)),
32 glm::dvec3(b_m * glm::dvec4(b_box.max.x, b_box.min.y, b_box.max.z, 1.0)),
33 glm::dvec3(b_m * glm::dvec4(b_box.max.x, b_box.max.y, b_box.min.z, 1.0)),
34 glm::dvec3(b_m * glm::dvec4(b_box.max.x, b_box.max.y, b_box.max.z, 1.0)),
37 glm::dvec3 axes[15] = {
44 glm::normalize(glm::cross(glm::dvec3(a_m[0]), glm::dvec3(b_m[0]))),
45 glm::normalize(glm::cross(glm::dvec3(a_m[0]), glm::dvec3(b_m[1]))),
46 glm::normalize(glm::cross(glm::dvec3(a_m[0]), glm::dvec3(b_m[2]))),
47 glm::normalize(glm::cross(glm::dvec3(a_m[1]), glm::dvec3(b_m[0]))),
48 glm::normalize(glm::cross(glm::dvec3(a_m[1]), glm::dvec3(b_m[1]))),
49 glm::normalize(glm::cross(glm::dvec3(a_m[1]), glm::dvec3(b_m[2]))),
50 glm::normalize(glm::cross(glm::dvec3(a_m[2]), glm::dvec3(b_m[0]))),
51 glm::normalize(glm::cross(glm::dvec3(a_m[2]), glm::dvec3(b_m[1]))),
52 glm::normalize(glm::cross(glm::dvec3(a_m[2]), glm::dvec3(b_m[2]))),
55 depth = std::numeric_limits<double>::infinity();
59 for (const glm::dvec3 &axis : axes) {
60 if (glm::any(glm::isnan(axis))) {
61 // can result from the cross products if A and B have parallel axes
65 double a_min = std::numeric_limits<double>::infinity();
66 double a_max = -std::numeric_limits<double>::infinity();
67 for (const glm::dvec3 &corner : a_corners) {
68 double val = glm::dot(corner, axis);
69 a_min = std::min(a_min, val);
70 a_max = std::max(a_max, val);
73 double b_min = std::numeric_limits<double>::infinity();
74 double b_max = -std::numeric_limits<double>::infinity();
75 for (const glm::dvec3 &corner : b_corners) {
76 double val = glm::dot(corner, axis);
77 b_min = std::min(b_min, val);
78 b_max = std::max(b_max, val);
81 if (a_max < b_min || b_max < a_min) return false;
83 double overlap = std::min(a_max, b_max) - std::max(a_min, b_min);
84 if (overlap < depth) {
92 normal = axes[min_axis];
105 double t_max = std::numeric_limits<double>::infinity();
106 const glm::dvec3 aabb_pos(M[3].x, M[3].y, M[3].z);
107 const glm::dvec3 delta = aabb_pos - ray.Origin();
109 glm::dvec3 t1(t_min, t_min, t_min), t2(t_max, t_max, t_max);
111 for (int i = 0; i < 3; ++i) {
112 const glm::dvec3 axis(M[i].x, M[i].y, M[i].z);
113 const double e = glm::dot(axis, delta);
114 const double f = glm::dot(axis, ray.Direction());
116 if (std::abs(f) > std::numeric_limits<double>::epsilon()) {
117 t1[i] = (e + aabb.min[i]) / f;
118 t2[i] = (e + aabb.max[i]) / f;
120 t_min = std::max(t_min, std::min(t1[i], t2[i]));
121 t_max = std::min(t_max, std::max(t1[i], t2[i]));
127 if (aabb.min[i] - e > 0.0 || aabb.max[i] - e < 0.0) {
135 glm::dvec3 min_all(glm::min(t1, t2));
136 if (min_all.x > min_all.y) {
137 if (min_all.x > min_all.z) {
138 normal = glm::dvec3(t2.x < t1.x ? 1 : -1, 0, 0);
140 normal = glm::dvec3(0, 0, t2.z < t1.z ? 1 : -1);
142 } else if (min_all.y > min_all.z) {
143 normal = glm::dvec3(0, t2.y < t1.y ? 1 : -1, 0);
145 normal = glm::dvec3(0, 0, t2.z < t1.z ? 1 : -1);
157 const glm::dvec3 diff(s.origin - r.Origin());
158 if (glm::dot(diff, r.Direction()) < 0.0) {
159 if (glm::length2(diff) > s.radius * s.radius) return false;
160 if (std::abs(glm::length2(diff) - s.radius * s.radius) < std::numeric_limits<double>::epsilon() * s.radius) {
161 normal = glm::normalize(-diff);
165 const glm::dvec3 pc(r.Direction() * glm::dot(r.Direction(), diff) + r.Origin());
166 double idist = std::sqrt(s.radius * s.radius - glm::length2(pc - s.origin));
167 dist = idist - glm::length(pc - r.Origin());
168 normal = glm::normalize((r.Origin() + (r.Direction() * dist)) - s.origin);
171 const glm::dvec3 pc(r.Direction() * glm::dot(r.Direction(), diff) + r.Origin());
172 if (glm::length2(s.origin - pc) > s.radius * s.radius) return false;
173 double idist = std::sqrt(s.radius * s.radius - glm::length2(pc - s.origin));
174 if (glm::length2(diff) > s.radius * s.radius) {
175 dist = glm::length(pc - r.Origin()) - idist;
177 dist = glm::length(pc - r.Origin()) + idist;
179 normal = glm::normalize((r.Origin() + (r.Direction() * dist)) - s.origin);