cgv
Loading...
Searching...
No Matches
intersection.h
1#pragma once
2
3#include "functions.h"
4#include "fvec.h"
5#include "pose.h"
6#include "ray.h"
7#include <limits>
8
10namespace cgv {
11namespace math {
12
22template <typename T>
23int ray_box_intersection(const ray<T, 3>& ray, fvec<T, 3> extent, fvec<T, 2>& out_ts, fvec<T, 3>* out_normal = nullptr) {
24
25 fvec<T, 3> m = fvec<T, 3>(T(1)) / ray.direction; // could be precomputed if traversing a set of aligned boxes
26 fvec<T, 3> n = m * ray.origin; // could be precomputed if traversing a set of aligned boxes
27 fvec<T, 3> k = abs(m) * extent;
28 fvec<T, 3> t1 = -n - k;
29 fvec<T, 3> t2 = -n + k;
30 T t_near = std::max(std::max(t1.x(), t1.y()), t1.z());
31 T t_far = std::min(std::min(t2.x(), t2.y()), t2.z());
32
33 if(t_near > t_far || t_far < T(0))
34 return 0;
35
36 out_ts[0] = t_near;
37 out_ts[1] = t_far;
38
39 if(out_normal)
40 *out_normal = -sign(ray.direction)
41 * step(fvec<T, 3>(t1.y(), t1.z(), t1.x()), fvec<T, 3>(t1.x(), t1.y(), t1.z()))
42 * step(fvec<T, 3>(t1.z(), t1.x(), t1.y()), fvec<T, 3>(t1.x(), t1.y(), t1.z()));
43
44 return 2;
45}
46
56template <typename T>
57int ray_box_intersection(const ray<T, 3> &ray, const fvec<T, 3> &min, const fvec<T, 3> &max, fvec<T, 2>& out_ts) {
58
59 fvec<T, 3> t0 = (min - ray.origin) / ray.direction;
60 fvec<T, 3> t1 = (max - ray.origin) / ray.direction;
61
62 if(t0.x() > t1.x())
63 std::swap(t0.x(), t1.x());
64
65 if(t0.y() > t1.y())
66 std::swap(t0.y(), t1.y());
67
68 if(t0.z() > t1.z())
69 std::swap(t0.z(), t1.z());
70
71 if(t0.x() > t1.y() || t0.y() > t1.x() ||
72 t0.x() > t1.z() || t0.z() > t1.x() ||
73 t0.z() > t1.y() || t0.y() > t1.z())
74 return 0;
75
76 T t_near = std::max(std::max(t0.x(), t0.y()), t0.z());
77 T t_far = std::min(std::min(t1.x(), t1.y()), t1.z());
78
79 if(t_near > t_far)
80 std::swap(t_near, t_far);
81
82 out_ts[0] = t_near;
83 out_ts[1] = t_far;
84
85 return 2;
86}
87
99template <typename T>
100int ray_cylinder_intersection(const ray<T, 3>& ray, const fvec<T, 3>& position, const fvec<T, 3>& axis, T radius, T& out_t, fvec<T, 3>* out_normal = nullptr) {
101
102 fvec<T, 3> oc = ray.origin - position;
103 T caca = dot(axis, axis);
104 T card = dot(axis, ray.direction);
105 T caoc = dot(axis, oc);
106 T a = caca - card * card;
107 T b = caca * dot(oc, ray.direction) - caoc * card;
108 T c = caca * dot(oc, oc) - caoc * caoc - radius * radius * caca;
109 T h = b * b - a * c;
110
111 if(h < T(0))
112 return 0;
113
114 h = std::sqrt(h);
115 out_t = (-b - h) / a;
116
117 // body
118 T y = caoc + out_t * card;
119 if(y > T(0) && y < caca) {
120 if(out_normal)
121 *out_normal = (oc + out_t * ray.direction - axis * y / caca) / radius;
122 return 1;
123 }
124
125 // caps
126 out_t = ((y < T(0) ? T(0) : caca) - caoc) / card;
127 if(std::abs(b + a * out_t) < h) {
128 if(out_normal)
129 *out_normal = axis * sign(y) / caca;
130 return 1;
131 }
132
133 return 0;
134}
135
145template <typename T>
146int ray_cylinder_intersection2(const ray<T, 3>& ray, const fvec<T, 3>& start_position, const fvec<T, 3>& end_position, T radius, T& out_t, fvec<T, 3>* out_normal = nullptr) {
147
148 return ray_cylinder_intersection(ray, start_position, end_position - start_position, radius, out_t, out_normal);
149}
150
158template <typename T>
159int ray_plane_intersection(const ray<T, 3>& ray, const fvec<T, 3>& origin, const fvec<T, 3>& normal, T& out_t) {
160
161 T denom = dot(normal, ray.direction);
162 if(std::abs(denom) < std::numeric_limits<T>::epsilon())
163 return 0;
164
165 out_t = dot(origin - ray.origin, normal) / denom;
166 return 1;
167};
168
178template <typename T>
179int ray_axis_aligned_rectangle_intersection(const ray<T, 3>& ray, const fvec<T, 3>& position, const fvec<T, 2>& extent, int axis_index, T& out_t, fvec<T, 2>* out_uv = nullptr) {
180
181 assert(axis_index >= 0 && axis_index < 3);
182
183 fvec<T, 3> normal = { T(0) };
184 normal[axis_index] = T(1);
185
186 T t = std::numeric_limits<T>::max();
187 if(cgv::math::ray_plane_intersection(ray, position, normal, t)) {
188 fvec<T, 3> intersection_position = ray.position(t);
189 intersection_position -= position;
190
191 vec2 uv;
192 switch(axis_index) {
193 case 0:
194 uv[0] = intersection_position[1];
195 uv[1] = intersection_position[2];
196 break;
197 case 1:
198 uv[0] = intersection_position[0];
199 uv[1] = intersection_position[2];
200 break;
201 case 2:
202 uv[0] = intersection_position[0];
203 uv[1] = intersection_position[1];
204 break;
205 default:
206 return 0;
207 }
208
209 uv += T(0.5) * extent;
210
211 if(uv[0] >= T(0) && uv[0] <= extent.x() && uv[1] >= T(0) && uv[1] <= extent.y()) {
212 out_t = t;
213 if(out_uv)
214 *out_uv = uv / extent;
215 return 1;
216 }
217 }
218
219 return 0;
220}
221
232template <typename T>
233int ray_parallelogram_intersection(const ray<T, 3>& ray, const fvec<T, 3>& origin, const fvec<T, 3> edge_u, const fvec<T, 3>& edge_v, T& out_t, fvec<T, 3>* out_normal = nullptr, fvec<T, 2>* out_uv = nullptr) {
234
235 fvec<T, 3> normal = normalize(cross(edge_u, edge_v));
236
237 T sf = T(0);
238 int ku = 0;
239 int kv = 1;
240
241 // decide on best projection plane based on projected surface area
242 //area in xy plane
243 T axy = edge_u.x() * edge_u.x() + edge_u.y() * edge_u.y();
244 axy *= edge_v.x() * edge_v.x() + edge_v.y() * edge_v.y();
245
246 //area in xz plane
247 T axz = edge_u.x() * edge_u.x() + edge_u.z() * edge_u.z();
248 axz *= edge_v.x() * edge_v.x() + edge_v.z() * edge_v.z();
249
250 //area in yz plane
251 T ayz = edge_u.y() * edge_u.y() + edge_u.z() * edge_u.z();
252 ayz *= edge_v.y() * edge_v.y() + edge_v.z() * edge_v.z();
253
254 if(axy > axz) {
255 if(axy > ayz) {
256 //xy
257 ku = 0;
258 kv = 1;
259 sf = normal.z() < T(0) ? T(1) : -T(1);
260 } else {
261 //yz
262 ku = 1;
263 kv = 2;
264 sf = normal.x() < T(0) ? T(1) : -T(1);
265 }
266 } else {
267 if(axz > ayz) {
268 //xz
269 ku = 0;
270 kv = 2;
271 sf = normal.y() < T(0) ? -T(1) : T(1);
272 } else {
273 //yz
274 ku = 1;
275 kv = 2;
276 sf = normal.x() < T(0) ? T(1) : -T(1);
277 }
278 }
279
280 T ndd = dot(normal, ray.direction);
281 if(std::abs(ndd) < std::numeric_limits<T>::epsilon())
282 return 0;
283
284 T t = dot(normal, origin - ray.origin) / ndd;
285
286 //ray intersects plane
287 //now test if hitpoint is inside parallelogram
288 fvec<T, 3> x = ray.position(t);
289 fvec<T, 2> x2d(x[ku] - origin[ku], x[kv] - origin[kv]);
290
291 fvec<T, 2> e1(edge_u[ku], edge_u[kv]);
292 fvec<T, 2> e2(edge_v[ku], edge_v[kv]);
293
294 T s = e1.x() * x2d.y() - e1.y() * x2d.x();
295 if(sf * s > -std::numeric_limits<T>::epsilon())
296 return 0;
297
298 s = e2.x() * x2d.y() - e2.y() * x2d.x();
299 if(sf * s < std::numeric_limits<T>::epsilon())
300 return 0;
301
302 x2d -= (e1 + e2);
303
304 s = e1.y() * x2d.x() - e1.x() * x2d.y();
305 if(sf * s > -std::numeric_limits<T>::epsilon())
306 return 0;
307
308 s = e2.y() * x2d.x() - e2.x() * x2d.y();
309 if(sf * s < std::numeric_limits<T>::epsilon())
310 return 0;
311
312 out_t = t;
313
314 if(out_normal)
315 *out_normal = normal;
316
317 if(out_uv) {
318 fvec<T, 2> uv = x2d;
319 uv.x() /= length(e1);
320 uv.y() /= length(e2);
321 *out_uv = uv;
322 }
323
324 return 1;
325};
326
337template <typename T>
338int ray_rectangle_intersection(const ray<T, 3>& ray, const fvec<T, 3>& position, const fvec<T, 2> extent, const quaternion<T>& rotation, T& out_t, fvec<T, 3>* out_normal = nullptr, fvec<T, 2>* out_uv = nullptr) {
339
340 // define tangent and bitangent assuming the normal is (0, 1, 0) without rotation
341 fvec<T, 3> tangent = { T(1), T(0), T(0) };
342 fvec<T, 3> bitangent = { T(0), T(1), T(0) };
343
344 tangent = rotation.apply(tangent);
345 bitangent = rotation.apply(bitangent);
346
347 fvec<T, 3> corner = position - T(0.5) * extent.x() * tangent - T(0.5) * extent.y() * bitangent;
348
349 fvec<T, 3> edge_u = extent.x() * tangent;
350 fvec<T, 3> edge_v = extent.y() * bitangent;
351
352 return ray_parallelogram_intersection(ray, corner, edge_u, edge_v, out_t, out_normal, out_uv);
353};
354
362template <typename T>
363int ray_sphere_intersection(const ray<T, 3>& ray, const fvec<T, 3>& center, T radius, fvec<T, 2>& out_ts) {
364
365 fvec<T, 3> d = ray.origin - center;
366 T il = T(1) / dot(ray.direction, ray.direction);
367 T b = il * dot(d, ray.direction);
368 T c = il * (dot(d, d) - radius * radius);
369 T D = b * b - c;
370
371 if(D < T(0))
372 return 0;
373
374 if(D < std::numeric_limits<T>::epsilon()) {
375 out_ts = -b;
376 return 1;
377 }
378
379 D = std::sqrt(D);
380 out_ts[0] = -b - D;
381 out_ts[1] = -b + D;
382
383 return 2;
384}
385
394template <typename T>
395int first_ray_sphere_intersection(const ray<T, 3>& ray, const fvec<T, 3>& center, T radius, T& out_t, fvec<T, 3>* out_normal = nullptr) {
396
397 fvec<T, 2> ts;
398 int k = ray_sphere_intersection(ray, center, radius, ts);
399
400 if(k == 1 || (k == 2 && ts[0] > T(0)))
401 out_t = ts[0];
402 else if(k == 2 && ts[1] > T(0))
403 out_t = ts[1];
404 else
405 return 0;
406
407 if(out_normal)
408 *out_normal = normalize(ray.position(out_t) - center);
409
410 return 1;
411}
412
421template <typename T>
422int ray_torus_intersection(const ray<T, 3>& ray, T large_radius, T small_radius, T& out_t, fvec<T, 3>* out_normal = nullptr) {
423
424 T po = T(1);
425 T Ra2 = large_radius * large_radius;
426 T ra2 = small_radius * small_radius;
427 T m = dot(ray.origin, ray.origin);
428 T n = dot(ray.origin, ray.direction);
429 T k = (m + Ra2 - ra2) / T(2);
430 T k3 = n;
431 const fvec<T, 2>& ro_xy = reinterpret_cast<const fvec<T, 2>&>(ray.origin);
432 const fvec<T, 2>& rd_xy = reinterpret_cast<const fvec<T, 2>&>(ray.direction);
433 T k2 = n * n - Ra2 * dot(rd_xy, rd_xy) + k;
434 T k1 = n * k - Ra2 * dot(rd_xy, ro_xy);
435 T k0 = k * k - Ra2 * dot(ro_xy, ro_xy);
436
437 if(std::abs(k3 * (k3 * k3 - k2) + k1) < T(0.01)) {
438 po = T(-1);
439 T tmp = k1; k1 = k3; k3 = tmp;
440 k0 = T(1) / k0;
441 k1 = k1 * k0;
442 k2 = k2 * k0;
443 k3 = k3 * k0;
444 }
445
446 T c2 = k2 * T(2) - T(3) * k3 * k3;
447 T c1 = k3 * (k3 * k3 - k2) + k1;
448 T c0 = k3 * (k3 * (c2 + T(2) * k2) - T(8) * k1) + T(4) * k0;
449 c2 /= T(3);
450 c1 *= T(2);
451 c0 /= T(3);
452 T Q = c2 * c2 + c0;
453 T R = c2 * c2 * c2 - T(3) * c2 * c0 + c1 * c1;
454 T h = R * R - Q * Q * Q;
455
456 if(h >= T(0)) {
457 h = std::sqrt(h);
458 T v = sign(R + h) * std::pow(std::abs(R + h), T(1) / T(3)); // cube root
459 T u = sign(R - h) * std::pow(std::abs(R - h), T(1) / T(3)); // cube root
460 fvec<T, 2> s = fvec<T, 2>((v + u) + T(4) * c2, (v - u) * std::sqrt(T(3)));
461 T y = std::sqrt(T(0.5) * (length(s) + s.x()));
462 T x = T(0.5) * s.y() / y;
463 T r = T(2) * c1 / (x * x + y * y);
464 T t1 = x - r - k3; t1 = (po < T(0)) ? T(2) / t1 : t1;
465 T t2 = -x - r - k3; t2 = (po < T(0)) ? T(2) / t2 : t2;
466
467 if(t1 > T(0)) out_t = t1;
468 if(t2 > T(0)) out_t = std::min(out_t, t2);
469
470 if(out_normal) {
471 fvec<T, 3> pos = ray.position(out_t);
472 *out_normal = normalize(pos * ((dot(pos, pos) - ra2) * fvec<T, 3>(T(1)) - Ra2 * fvec<T, 3>(T(1), T(1), T(-1))));
473 }
474
475 return 1; // 2
476 }
477
478 T sQ = std::sqrt(Q);
479 T w = sQ * cos(acos(-R / (sQ * Q)) / T(3));
480 T d2 = -(w + c2);
481
482 if (d2 < T(0))
483 return 0;
484
485 T d1 = std::sqrt(d2);
486 T h1 = std::sqrt(w - T(2) * c2 + c1 / d1);
487 T h2 = std::sqrt(w - T(2) * c2 - c1 / d1);
488 T t1 = -d1 - h1 - k3; t1 = (po < T(0)) ? T(2) / t1 : t1;
489 T t2 = -d1 + h1 - k3; t2 = (po < T(0)) ? T(2) / t2 : t2;
490 T t3 = d1 - h2 - k3; t3 = (po < T(0)) ? T(2) / t3 : t3;
491 T t4 = d1 + h2 - k3; t4 = (po < T(0)) ? T(2) / t4 : t4;
492
493 if (t1 > T(0)) out_t = t1;
494 if (t2 > T(0)) out_t = std::min(out_t, t2);
495 if (t3 > T(0)) out_t = std::min(out_t, t3);
496 if (t4 > T(0)) out_t = std::min(out_t, t4);
497
498 if(out_normal) {
499 fvec<T, 3> pos = ray.position(out_t);
500 *out_normal = normalize(pos * ((dot(pos, pos) - ra2) * fvec<T, 3>(T(1)) - Ra2 * fvec<T, 3>(T(1), T(1), T(-1))));
501 }
502
503 return 1; // 4
504}
505
516template <typename T>
517int ray_torus_intersection(const ray<T, 3>& ray, const fvec<T, 3>& center, const fvec<T, 3>& normal, T large_radius, T small_radius, T& out_t, fvec<T, 3>* out_normal = nullptr) {
518
519 // compute pose transformation
520 fmat<T, 3, 4> pose;
521 cgv::math::pose_position(pose) = center;
522 fvec<T, 3>& x = reinterpret_cast<fvec<T, 3>&>(pose[0]);
523 fvec<T, 3>& y = reinterpret_cast<fvec<T, 3>&>(pose[3]);
524 fvec<T, 3>& z = reinterpret_cast<fvec<T, 3>&>(pose[6]);
525 z = normal;
526 x = normal;
527 int i = std::abs(normal[0]) < std::abs(normal[1]) ? 0 : 1;
528 i = std::abs(normal[i]) < std::abs(normal[2]) ? i : 2;
529 x[i] = T(1);
530 y = normalize(cross(normal, x));
531 x = cross(y, normal);
532
533 cgv::math::ray<T, 3> transformed_ray;
534 transformed_ray.origin = cgv::math::inverse_pose_transform_point(pose, ray.origin);
535 transformed_ray.direction = cgv::math::inverse_pose_transform_vector(pose, ray.direction);
536
537 // transform ray into torus pose
538 int res = ray_torus_intersection(transformed_ray, large_radius, small_radius, out_t, out_normal);
539
540 // in case of intersection, transform normal back to world space
541 if(res)
542 *out_normal = cgv::math::pose_transform_vector(pose, *out_normal);
543
544 return res;
545}
546
547} // namespace math
548} // namespace cgv
cgv::math::ray
This class defines a template for n-dimensional rays with arbitrary data type defined by origin and d...
Definition ray.h:14
cgv
the cgv namespace
Definition print.h:11
cgv::vec2
cgv::math::fvec< float, 2 > vec2
declare type of 2d single precision floating point vectors
Definition fvec.h:667
pose.h
helper functions to work with poses that can be represented with 3x4 matrix or quaternion plus vector
cgv::math::pose_transform_vector
fvec< T, 3 > pose_transform_vector(const fmat< T, 3, 4 > &pose, const fvec< T, 3 > &v)
transform vector with pose matrix
Definition pose.h:30
cgv::math::pose_position
fvec< T, 3 > & pose_position(fmat< T, 3, 4 > &pose)
extract position vector from pose matrix
Definition pose.h:19
cgv::math::inverse_pose_transform_vector
fvec< T, 3 > inverse_pose_transform_vector(const fmat< T, 3, 4 > &pose, const fvec< T, 3 > &v)
transform vector with inverse of pose matrix
Definition pose.h:36
cgv::math::inverse_pose_transform_point
fvec< T, 3 > inverse_pose_transform_point(const fmat< T, 3, 4 > &pose, const fvec< T, 3 > &p)
transform point with inverse of pose matrix
Definition pose.h:33