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
147template <typename T>
148int 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) {
149
150 return ray_cylinder_intersection(ray, start_position, end_position - start_position, radius, out_t, out_normal);
151}
152
162template <typename T>
163int ray_plane_intersection(const ray<T, 3>& ray, const fvec<T, 3>& origin, const fvec<T, 3>& normal, T& out_t) {
164
165 T denom = dot(normal, ray.direction);
166 if(std::abs(denom) < std::numeric_limits<T>::epsilon())
167 return 0;
168
169 out_t = dot(origin - ray.origin, normal) / denom;
170 return 1;
171};
172
184template <typename T>
185int 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) {
186
187 assert(axis_index >= 0 && axis_index < 3);
188
189 fvec<T, 3> normal = { T(0) };
190 normal[axis_index] = T(1);
191
192 T t = std::numeric_limits<T>::max();
193 if(cgv::math::ray_plane_intersection(ray, position, normal, t)) {
194 fvec<T, 3> intersection_position = ray.position(t);
195 intersection_position -= position;
196
197 vec2 uv;
198 switch(axis_index) {
199 case 0:
200 uv[0] = intersection_position[1];
201 uv[1] = intersection_position[2];
202 break;
203 case 1:
204 uv[0] = intersection_position[0];
205 uv[1] = intersection_position[2];
206 break;
207 case 2:
208 uv[0] = intersection_position[0];
209 uv[1] = intersection_position[1];
210 break;
211 default:
212 return 0;
213 }
214
215 uv += T(0.5) * extent;
216
217 if(uv[0] >= T(0) && uv[0] <= extent.x() && uv[1] >= T(0) && uv[1] <= extent.y()) {
218 out_t = t;
219 if(out_uv)
220 *out_uv = uv / extent;
221 return 1;
222 }
223 }
224
225 return 0;
226}
227
240template <typename T>
241int 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) {
242
243 fvec<T, 3> normal = normalize(cross(edge_u, edge_v));
244
245 T sf = T(0);
246 int ku = 0;
247 int kv = 1;
248
249 // decide on best projection plane based on projected surface area
250 //area in xy plane
251 T axy = edge_u.x() * edge_u.x() + edge_u.y() * edge_u.y();
252 axy *= edge_v.x() * edge_v.x() + edge_v.y() * edge_v.y();
253
254 //area in xz plane
255 T axz = edge_u.x() * edge_u.x() + edge_u.z() * edge_u.z();
256 axz *= edge_v.x() * edge_v.x() + edge_v.z() * edge_v.z();
257
258 //area in yz plane
259 T ayz = edge_u.y() * edge_u.y() + edge_u.z() * edge_u.z();
260 ayz *= edge_v.y() * edge_v.y() + edge_v.z() * edge_v.z();
261
262 if(axy > axz) {
263 if(axy > ayz) {
264 //xy
265 ku = 0;
266 kv = 1;
267 sf = normal.z() < T(0) ? T(1) : -T(1);
268 } else {
269 //yz
270 ku = 1;
271 kv = 2;
272 sf = normal.x() < T(0) ? T(1) : -T(1);
273 }
274 } else {
275 if(axz > ayz) {
276 //xz
277 ku = 0;
278 kv = 2;
279 sf = normal.y() < T(0) ? -T(1) : T(1);
280 } else {
281 //yz
282 ku = 1;
283 kv = 2;
284 sf = normal.x() < T(0) ? T(1) : -T(1);
285 }
286 }
287
288 T ndd = dot(normal, ray.direction);
289 if(std::abs(ndd) < std::numeric_limits<T>::epsilon())
290 return 0;
291
292 T t = dot(normal, origin - ray.origin) / ndd;
293
294 //ray intersects plane
295 //now test if hitpoint is inside parallelogram
296 fvec<T, 3> x = ray.position(t);
297 fvec<T, 2> x2d(x[ku] - origin[ku], x[kv] - origin[kv]);
298
299 fvec<T, 2> e1(edge_u[ku], edge_u[kv]);
300 fvec<T, 2> e2(edge_v[ku], edge_v[kv]);
301
302 T s = e1.x() * x2d.y() - e1.y() * x2d.x();
303 if(sf * s > -std::numeric_limits<T>::epsilon())
304 return 0;
305
306 s = e2.x() * x2d.y() - e2.y() * x2d.x();
307 if(sf * s < std::numeric_limits<T>::epsilon())
308 return 0;
309
310 x2d -= (e1 + e2);
311
312 s = e1.y() * x2d.x() - e1.x() * x2d.y();
313 if(sf * s > -std::numeric_limits<T>::epsilon())
314 return 0;
315
316 s = e2.y() * x2d.x() - e2.x() * x2d.y();
317 if(sf * s < std::numeric_limits<T>::epsilon())
318 return 0;
319
320 out_t = t;
321
322 if(out_normal)
323 *out_normal = normal;
324
325 if(out_uv) {
326 fvec<T, 2> uv = x2d;
327 uv.x() /= length(e1);
328 uv.y() /= length(e2);
329 *out_uv = uv;
330 }
331
332 return 1;
333};
334
347template <typename T>
348int 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) {
349
350 // define tangent and bitangent assuming the normal is (0, 1, 0) without rotation
351 fvec<T, 3> tangent = { T(1), T(0), T(0) };
352 fvec<T, 3> bitangent = { T(0), T(1), T(0) };
353
354 tangent = rotation.apply(tangent);
355 bitangent = rotation.apply(bitangent);
356
357 fvec<T, 3> corner = position - T(0.5) * extent.x() * tangent - T(0.5) * extent.y() * bitangent;
358
359 fvec<T, 3> edge_u = extent.x() * tangent;
360 fvec<T, 3> edge_v = extent.y() * bitangent;
361
362 return ray_parallelogram_intersection(ray, corner, edge_u, edge_v, out_t, out_normal, out_uv);
363};
364
374template <typename T>
375int ray_sphere_intersection(const ray<T, 3>& ray, const fvec<T, 3>& center, T radius, fvec<T, 2>& out_ts) {
376
377 fvec<T, 3> d = ray.origin - center;
378 T il = T(1) / dot(ray.direction, ray.direction);
379 T b = il * dot(d, ray.direction);
380 T c = il * (dot(d, d) - radius * radius);
381 T D = b * b - c;
382
383 if(D < T(0))
384 return 0;
385
386 if(D < std::numeric_limits<T>::epsilon()) {
387 out_ts = -b;
388 return 1;
389 }
390
391 D = std::sqrt(D);
392 out_ts[0] = -b - D;
393 out_ts[1] = -b + D;
394
395 return 2;
396}
397
408template <typename T>
409int 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) {
410
411 fvec<T, 2> ts;
412 int k = ray_sphere_intersection(ray, center, radius, ts);
413
414 if(k == 1 || (k == 2 && ts[0] > T(0)))
415 out_t = ts[0];
416 else if(k == 2 && ts[1] > T(0))
417 out_t = ts[1];
418 else
419 return 0;
420
421 if(out_normal)
422 *out_normal = normalize(ray.position(out_t) - center);
423
424 return 1;
425}
426
437template <typename T>
438int ray_torus_intersection(const ray<T, 3>& ray, T large_radius, T small_radius, T& out_t, fvec<T, 3>* out_normal = nullptr) {
439
440 T po = T(1);
441 T Ra2 = large_radius * large_radius;
442 T ra2 = small_radius * small_radius;
443 T m = dot(ray.origin, ray.origin);
444 T n = dot(ray.origin, ray.direction);
445 T k = (m + Ra2 - ra2) / T(2);
446 T k3 = n;
447 const fvec<T, 2>& ro_xy = reinterpret_cast<const fvec<T, 2>&>(ray.origin);
448 const fvec<T, 2>& rd_xy = reinterpret_cast<const fvec<T, 2>&>(ray.direction);
449 T k2 = n * n - Ra2 * dot(rd_xy, rd_xy) + k;
450 T k1 = n * k - Ra2 * dot(rd_xy, ro_xy);
451 T k0 = k * k - Ra2 * dot(ro_xy, ro_xy);
452
453 if(std::abs(k3 * (k3 * k3 - k2) + k1) < T(0.01)) {
454 po = T(-1);
455 T tmp = k1; k1 = k3; k3 = tmp;
456 k0 = T(1) / k0;
457 k1 = k1 * k0;
458 k2 = k2 * k0;
459 k3 = k3 * k0;
460 }
461
462 T c2 = k2 * T(2) - T(3) * k3 * k3;
463 T c1 = k3 * (k3 * k3 - k2) + k1;
464 T c0 = k3 * (k3 * (c2 + T(2) * k2) - T(8) * k1) + T(4) * k0;
465 c2 /= T(3);
466 c1 *= T(2);
467 c0 /= T(3);
468 T Q = c2 * c2 + c0;
469 T R = c2 * c2 * c2 - T(3) * c2 * c0 + c1 * c1;
470 T h = R * R - Q * Q * Q;
471
472 if(h >= T(0)) {
473 h = std::sqrt(h);
474 T v = sign(R + h) * std::pow(std::abs(R + h), T(1) / T(3)); // cube root
475 T u = sign(R - h) * std::pow(std::abs(R - h), T(1) / T(3)); // cube root
476 fvec<T, 2> s = fvec<T, 2>((v + u) + T(4) * c2, (v - u) * std::sqrt(T(3)));
477 T y = std::sqrt(T(0.5) * (length(s) + s.x()));
478 T x = T(0.5) * s.y() / y;
479 T r = T(2) * c1 / (x * x + y * y);
480 T t1 = x - r - k3; t1 = (po < T(0)) ? T(2) / t1 : t1;
481 T t2 = -x - r - k3; t2 = (po < T(0)) ? T(2) / t2 : t2;
482
483 if(t1 > T(0)) out_t = t1;
484 if(t2 > T(0)) out_t = std::min(out_t, t2);
485
486 if(out_normal) {
487 fvec<T, 3> pos = ray.position(out_t);
488 *out_normal = normalize(pos * ((dot(pos, pos) - ra2) * fvec<T, 3>(T(1)) - Ra2 * fvec<T, 3>(T(1), T(1), T(-1))));
489 }
490
491 return 1; // 2
492 }
493
494 T sQ = std::sqrt(Q);
495 T w = sQ * cos(acos(-R / (sQ * Q)) / T(3));
496 T d2 = -(w + c2);
497
498 if (d2 < T(0))
499 return 0;
500
501 T d1 = std::sqrt(d2);
502 T h1 = std::sqrt(w - T(2) * c2 + c1 / d1);
503 T h2 = std::sqrt(w - T(2) * c2 - c1 / d1);
504 T t1 = -d1 - h1 - k3; t1 = (po < T(0)) ? T(2) / t1 : t1;
505 T t2 = -d1 + h1 - k3; t2 = (po < T(0)) ? T(2) / t2 : t2;
506 T t3 = d1 - h2 - k3; t3 = (po < T(0)) ? T(2) / t3 : t3;
507 T t4 = d1 + h2 - k3; t4 = (po < T(0)) ? T(2) / t4 : t4;
508
509 if (t1 > T(0)) out_t = t1;
510 if (t2 > T(0)) out_t = std::min(out_t, t2);
511 if (t3 > T(0)) out_t = std::min(out_t, t3);
512 if (t4 > T(0)) out_t = std::min(out_t, t4);
513
514 if(out_normal) {
515 fvec<T, 3> pos = ray.position(out_t);
516 *out_normal = normalize(pos * ((dot(pos, pos) - ra2) * fvec<T, 3>(T(1)) - Ra2 * fvec<T, 3>(T(1), T(1), T(-1))));
517 }
518
519 return 1; // 4
520}
521
534template <typename T>
535int 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) {
536
537 // compute pose transformation
538 fmat<T, 3, 4> pose;
539 cgv::math::pose_position(pose) = center;
540 fvec<T, 3>& x = reinterpret_cast<fvec<T, 3>&>(pose[0]);
541 fvec<T, 3>& y = reinterpret_cast<fvec<T, 3>&>(pose[3]);
542 fvec<T, 3>& z = reinterpret_cast<fvec<T, 3>&>(pose[6]);
543 z = normal;
544 x = normal;
545 int i = std::abs(normal[0]) < std::abs(normal[1]) ? 0 : 1;
546 i = std::abs(normal[i]) < std::abs(normal[2]) ? i : 2;
547 x[i] = T(1);
548 y = normalize(cross(normal, x));
549 x = cross(y, normal);
550
551 cgv::math::ray<T, 3> transformed_ray;
552 transformed_ray.origin = cgv::math::inverse_pose_transform_point(pose, ray.origin);
553 transformed_ray.direction = cgv::math::inverse_pose_transform_vector(pose, ray.direction);
554
555 // transform ray into torus pose
556 int res = ray_torus_intersection(transformed_ray, large_radius, small_radius, out_t, out_normal);
557
558 // in case of intersection, transform normal back to world space
559 if(res)
560 *out_normal = cgv::math::pose_transform_vector(pose, *out_normal);
561
562 return res;
563}
564
565} // namespace math
566} // 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:668
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