intersection.hh 49.4 KB
Newer Older
1
2
#pragma once

3
#include <typed-geometry/detail/optional.hh>
Philip Trettner's avatar
Philip Trettner committed
4
5
#include <typed-geometry/feature/assert.hh>
#include <typed-geometry/functions/basic/scalar_math.hh>
6

7
#include <typed-geometry/types/objects/aabb.hh>
8
#include <typed-geometry/types/objects/box.hh>
9
#include <typed-geometry/types/objects/capsule.hh>
Philip Trettner's avatar
Philip Trettner committed
10
#include <typed-geometry/types/objects/cylinder.hh>
Aaron Grabowy's avatar
Aaron Grabowy committed
11
#include <typed-geometry/types/objects/ellipse.hh>
12
#include <typed-geometry/types/objects/halfspace.hh>
13
#include <typed-geometry/types/objects/hemisphere.hh>
Aaron Grabowy's avatar
Aaron Grabowy committed
14
15
#include <typed-geometry/types/objects/inf_cone.hh>
#include <typed-geometry/types/objects/inf_cylinder.hh>
16
#include <typed-geometry/types/objects/line.hh>
Philip Trettner's avatar
Philip Trettner committed
17
#include <typed-geometry/types/objects/plane.hh>
18
#include <typed-geometry/types/objects/pyramid.hh>
Philip Trettner's avatar
Philip Trettner committed
19
#include <typed-geometry/types/objects/ray.hh>
20
#include <typed-geometry/types/objects/segment.hh>
Philip Trettner's avatar
Philip Trettner committed
21
22
#include <typed-geometry/types/objects/sphere.hh>
#include <typed-geometry/types/objects/triangle.hh>
23

Philip Trettner's avatar
Philip Trettner committed
24
25
26
27
#include <typed-geometry/functions/vector/cross.hh>
#include <typed-geometry/functions/vector/dot.hh>
#include <typed-geometry/functions/vector/length.hh>

28
#include "aabb.hh"
29
#include "contains.hh"
30
#include "coordinates.hh"
Philip Trettner's avatar
Philip Trettner committed
31
#include "direction.hh"
32
#include "faces.hh"
33
#include "normal.hh"
34

Aaron Grabowy's avatar
Aaron Grabowy committed
35
36
#include <utility>

37
38
// family of intersection functions:

Philip Trettner's avatar
Philip Trettner committed
39
40
41
// intersects(a, b)              -> bool
// intersection(a, b)            -> ???
// intersection_safe(a, b)       -> optional<???>
42
// intersection_parameter(a, b)  -> coords? (for a line or a ray: hits<N, ScalarT> or optional<hit_interval> (when b is solid))
Philip Trettner's avatar
Philip Trettner committed
43
44
// intersection_parameters(a, b) -> pair<coords, coords>?
// intersection_exact(a, b)      -> variant
45
46
// closest_intersection(a, b)            -> position/object (for a ray: optional<pos>)
// closest_intersection_parameter(a, b)  -> coords (for a ray: optional<ScalarT>)
Philip Trettner's avatar
Philip Trettner committed
47
48
49
50

// "intersects" returns true iff any point lies in a and in b
// "intersection" returns an object describing the intersection (NOTE: does NOT handle degenerate cases)
// "intersection_safe" is the same as "intersection" but returns nullopt for degenerate cases
51
// "intersection_parameter" returns coordinates for the first object such that a[coords] == intersection(a, b)
Philip Trettner's avatar
Philip Trettner committed
52
53
// "intersection_parameters" returns coordinates for both objects
// "intersection_exact" returns a variant type describing all possible intersections, including degenerate cases
54
// the "closest_" variants only return the closest intersection for objects where that concept is applicable (e.g. for rays)
Philip Trettner's avatar
Philip Trettner committed
55
56
57
58
59

// Notes:
//  - intersection_exact is currently unsupported
//  - intersection_safe is currently unsupported
//  - for more elaborate ray-tracing, a future ray_cast function will exist (which also returns the intersection normal)
60

61
// Implementation guidelines:
62
63
64
65
66
67
68
69
70
71
72
// if object has boundary_of(obj) defined
//      explicit intersection_parameter(line, obj_boundary), which gives intersection_parameter(line, obj)
//      if obj is infinite and contains(obj, pos) is not cheap, intersection_parameter(line, obj) can be implemented additionally
// else
//      explicit intersection_parameter(line, obj)
// this gives intersection, intersects, closest_intersection_parameter, and closest_intersection for line and ray
//
// if closest ray intersection is faster than computing all line intersections
//      explicit closest_intersection_parameter(ray, obj), same for obj_boundary
// this is then also used by closest_intersection and intersects
//
73
// explicit intersects(obj, aabb), which gives intersects(aabb, obj)
74
75
//
// for convex compound objects (like cylinder or pyramids), decompose the object into primitive shapes and pass them to a helper function:
76
77
// - call merge_hits(line, objPart1, objPart2, ...) in the implementation of intersection_parameter(line, obj_boundary)
// - call intersects_any(lineRayObj, objPart1, objPart2, ...) in the implementation of intersects(lineRayObj, obj<TraitsT>), which can shortcut and be faster than the default
78

79

80
81
namespace tg
{
Philip Trettner's avatar
Philip Trettner committed
82
83
84
85
86
87
88
89
// ====================================== Result Structs ======================================

/// ordered list of ray intersection hits
/// behaves like a container with
///   .size()
///   operator[]
///   range-based-for
template <int MaxHits, class HitT>
90
struct hits
91
{
92
    static constexpr bool is_hits = true; // tag
93
    static constexpr int max_hits = MaxHits;
94
    using hit_t = HitT;
Philip Trettner's avatar
Philip Trettner committed
95
96

    template <class OtherT>
97
    using as_hits = hits<MaxHits, OtherT>;
Philip Trettner's avatar
Philip Trettner committed
98

99
100
    [[nodiscard]] int size() const { return _size; }
    [[nodiscard]] bool any() const { return _size > 0; }
Philip Trettner's avatar
Philip Trettner committed
101

102
    HitT const& operator[](int idx) const
Philip Trettner's avatar
Philip Trettner committed
103
104
105
106
    {
        TG_ASSERT(0 <= idx && idx < _size);
        return _hit[idx];
    }
107
108
109
110
111
    [[nodiscard]] HitT const& first() const
    {
        TG_ASSERT(_size > 0);
        return _hit[0];
    }
112
113
114
115
116
    [[nodiscard]] HitT const& last() const
    {
        TG_ASSERT(_size > 0);
        return _hit[_size - 1];
    }
Philip Trettner's avatar
Philip Trettner committed
117

118
119
    [[nodiscard]] HitT const* begin() const { return _hit; }
    [[nodiscard]] HitT const* end() const { return _hit + _size; }
Philip Trettner's avatar
Philip Trettner committed
120

121
122
    hits() = default;
    hits(HitT* hits, int size) : _size(size)
Philip Trettner's avatar
Philip Trettner committed
123
124
125
126
    {
        for (auto i = 0; i < size; ++i)
            _hit[i] = hits[i];
    }
127
    template <typename... HitTs>
128
    hits(HitTs... hits) : _size(sizeof...(HitTs)), _hit{hits...}
129
130
    {
    }
Philip Trettner's avatar
Philip Trettner committed
131
132
133

private:
    int _size = 0;
Aaron Grabowy's avatar
Aaron Grabowy committed
134
    HitT _hit[MaxHits] = {};
135
136
};

137
/// describes a continuous interval on a line or ray between start and end
138
template <class ScalarT>
139
struct hit_interval
140
{
141
    static constexpr bool is_hit_interval = true; // tag
142

143
144
145
    ScalarT start;
    ScalarT end;

146
    [[nodiscard]] constexpr bool is_unbounded() const { return end == tg::max<ScalarT>() || start == tg::min<ScalarT>(); }
147
148
};

Philip Trettner's avatar
Philip Trettner committed
149

150
151
152
153
154
155
// ====================================== Helper functions ======================================

namespace detail
{
// intersects the given line with all given objects and returns the concatenated intersections. A maximal number of 2 intersections is assumed.
template <int D, class ScalarT, class... Objs>
156
[[nodiscard]] constexpr hits<2, ScalarT> merge_hits(line<D, ScalarT> const& line, Objs const&... objs)
157
{
Aaron Grabowy's avatar
Aaron Grabowy committed
158
    ScalarT hits[2] = {};
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
    hits[0] = tg::max<ScalarT>();
    hits[1] = tg::min<ScalarT>();
    auto numHits = 0;

    const auto find_hits = [&](const auto& obj) {
        const auto inters = intersection_parameter(line, obj);
        for (const auto& inter : inters)
        {
            hits[0] = tg::min(hits[0], inter);
            hits[1] = tg::max(hits[1], inter);
            numHits++;
        }
    };
    (find_hits(objs), ...);

    TG_ASSERT(numHits <= 2);
    return {hits, numHits};
}
Aaron Grabowy's avatar
Aaron Grabowy committed
177
178
template <int D, class ScalarT, class ObjT, u64 N, size_t... I>
[[nodiscard]] constexpr hits<2, ScalarT> merge_hits_array(line<D, ScalarT> const& line, array<ObjT, N> objs, std::index_sequence<I...>)
179
180
181
{
    return merge_hits(line, objs[I]...);
}
182
183
184
185
186
187
188
189

// returns true, iff the given line or ray object intersects any of the given other objects (with short-circuiting after the first intersection)
template <class Obj, class... Objs>
[[nodiscard]] constexpr bool intersects_any(Obj const& obj, Objs const&... objs)
{
    return (intersects(obj, objs) || ...);
}

Aaron Grabowy's avatar
Aaron Grabowy committed
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
// Solves the quadratic equation ax^2 + bx + c = 0
template <class ScalarT>
[[nodiscard]] constexpr hits<2, ScalarT> solve_quadratic(ScalarT const& a, ScalarT const& b, ScalarT const& c)
{
    const auto discriminant = b * b - ScalarT(4) * a * c;
    if (discriminant < ScalarT(0))
        return {}; // No solution

    const auto sqrtD = sqrt(discriminant);
    const auto t1 = (-b - sqrtD) / (ScalarT(2) * a);
    const auto t2 = (-b + sqrtD) / (ScalarT(2) * a);

    const auto [tMin, tMax] = minmax(t1, t2);
    return {tMin, tMax};
}

206
template <class A, class B>
207
using try_closest_intersection_parameter = decltype(closest_intersection_parameter(std::declval<A const&>(), std::declval<B const&>()));
208
209
210
}


Philip Trettner's avatar
Philip Trettner committed
211
// ====================================== Default Implementations ======================================
Philip Trettner's avatar
Philip Trettner committed
212
// TODO: intersection_parameter from intersection_parameters
Philip Trettner's avatar
Philip Trettner committed
213
214

// returns whether two objects intersect
215
// if intersection is available and applicable, use that
Philip Trettner's avatar
Philip Trettner committed
216
template <class A, class B>
217
218
[[nodiscard]] constexpr auto intersects(A const& a, B const& b)
    -> std::enable_if_t<!can_apply<detail::try_closest_intersection_parameter, A, B>, decltype(intersection(a, b).has_value())>
219
{
Philip Trettner's avatar
Philip Trettner committed
220
221
    return intersection(a, b).has_value();
}
222

223
// if closest intersection parameter is available and applicable, use that
224
225
template <class A, class B>
[[nodiscard]] constexpr auto intersects(A const& a, B const& b) -> decltype(closest_intersection_parameter(a, b).has_value())
226
{
227
    return closest_intersection_parameter(a, b).has_value();
228
229
}

230
231
232
233
234
235
236
237
238
239
240
241
242
// parameters for intersects with aabb can switch order
template <int D, class ScalarT, class Obj>
[[nodiscard]] constexpr bool intersects(aabb<D, ScalarT> const& b, Obj const& obj)
{
    return intersects(obj, b);
}
// Explicit intersects aabb aabb to prevent infinite recursion
template <int D, class ScalarT>
[[nodiscard]] constexpr bool intersects(aabb<D, ScalarT> const& a, aabb<D, ScalarT> const& b)
{
    return intersection(a, b).has_value();
}

Philip Trettner's avatar
Philip Trettner committed
243
244
// if a value-typed intersection parameter is available and applicable, use that
template <class A, class B>
245
[[nodiscard]] constexpr auto intersection(A const& a, B const& b) -> decltype(a[intersection_parameter(a, b)])
246
{
Philip Trettner's avatar
Philip Trettner committed
247
248
    return a[intersection_parameter(a, b)];
}
249

Philip Trettner's avatar
Philip Trettner committed
250
251
// if an optional intersection parameter is available and applicable, use that
template <class A, class B>
252
[[nodiscard]] constexpr auto intersection(A const& a, B const& b) -> optional<decltype(a[intersection_parameter(a, b).value()])>
253
{
Philip Trettner's avatar
Philip Trettner committed
254
255
256
257
    if (auto t = intersection_parameter(a, b); t.has_value())
        return a[t.value()];
    return {};
}
258

Philip Trettner's avatar
Philip Trettner committed
259
260
// if a value-typed closest intersection parameter is available and applicable, use that
template <class A, class B>
261
[[nodiscard]] constexpr auto closest_intersection(A const& a, B const& b) -> decltype(a[closest_intersection_parameter(a, b)])
262
{
Philip Trettner's avatar
Philip Trettner committed
263
264
    return a[closest_intersection_parameter(a, b)];
}
265

Philip Trettner's avatar
Philip Trettner committed
266
// if an optional closest intersection parameter is available and applicable, use that
267
template <class A, class B>
268
[[nodiscard]] constexpr auto closest_intersection(A const& a, B const& b) -> optional<decltype(a[closest_intersection_parameter(a, b).value()])>
269
{
Philip Trettner's avatar
Philip Trettner committed
270
271
272
    if (auto t = closest_intersection_parameter(a, b); t.has_value())
        return a[t.value()];
    return {};
273
274
}

275
// if hits intersection parameter is available, use that
276
277
template <class A, class B>
[[nodiscard]] constexpr auto intersection(A const& a, B const& b) ->
278
    typename decltype(intersection_parameter(a, b))::template as_hits<typename A::pos_t>
Philip Trettner's avatar
Philip Trettner committed
279
{
280
    auto ts = intersection_parameter(a, b);
Aaron Grabowy's avatar
Aaron Grabowy committed
281
    typename A::pos_t hits[ts.max_hits] = {};
Philip Trettner's avatar
Philip Trettner committed
282
    for (auto i = 0; i < ts.size(); ++i)
283
        hits[i] = a[ts[i]];
Philip Trettner's avatar
Philip Trettner committed
284
285
286
    return {hits, ts.size()};
}

287
// if an optional hit_interval intersection parameter is available, use that
288
289
template <class A, class B>
[[nodiscard]] constexpr auto intersection(A const& a, B const& b)
290
    -> std::enable_if<decltype(intersection_parameter(a, b).value())::is_hit_interval, // condition
291
    optional<segment<object_traits<A>::domain_dimension, typename object_traits<A>::scalar_t>>> // return type
292
{
293
    auto ts = intersection_parameter(a, b);
294
    if (ts.has_value())
295
        return {{a[ts.value().start], a[ts.value().end]}};
296
297
298
    return {};
}

299
// if hits intersection parameter is available, use that
300
301
template <class A, class B>
[[nodiscard]] constexpr auto closest_intersection_parameter(A const& a, B const& b) -> optional<typename decltype(intersection_parameter(a, b))::hit_t>
Philip Trettner's avatar
Philip Trettner committed
302
{
303
    const auto hits = intersection_parameter(a, b);
Philip Trettner's avatar
Philip Trettner committed
304
    if (hits.any())
305
        return hits.first();
Philip Trettner's avatar
Philip Trettner committed
306
307
308
    return {};
}

309
// if an optional hit_interval intersection parameter is available, use that
310
311
template <class A, class B>
[[nodiscard]] constexpr auto closest_intersection_parameter(A const& a, B const& b) -> optional<decltype(intersection_parameter(a, b).value().start)>
312
{
313
    const auto hits = intersection_parameter(a, b);
314
315
316
317
318
    if (hits.has_value())
        return hits.value().start;
    return {};
}

319
// if boundary_of a solid object returns hits, use this to construct the hit_interval result of the solid intersection
320
template <int D, class ScalarT, class Obj>
321
[[nodiscard]] constexpr auto intersection_parameter(line<D, ScalarT> const& l, Obj const& obj)
322
    -> enable_if<!std::is_same_v<Obj, decltype(boundary_of(obj))>, optional<hit_interval<ScalarT>>>
323
{
324
    const hits<2, ScalarT> inter = intersection_parameter(l, boundary_of(obj));
325

326
327
    if (inter.size() == 2)
        return {{inter[0], inter[1]}};
328

329
    if constexpr (object_traits<Obj>::is_finite)
330
    {
331
332
        TG_ASSERT(inter.size() == 0); // a line intersects a finite solid object either twice or not at all
        return {};
333
    }
334
335
336
337
338
339
340
341
    else
    {
        if (inter.size() == 0)
        {
            if (contains(obj, l.pos))
                return {{tg::min<ScalarT>(), tg::max<ScalarT>()}};
            return {};
        }
342

343
344
345
346
347
        TG_ASSERT(inter.size() == 1); // no other number remains
        if (contains(obj, l[inter.first() + ScalarT(1)])) // the line continues into the object
            return {{inter.first(), tg::max<ScalarT>()}};
        return {{tg::min<ScalarT>(), inter.first()}};
    }
348
349
}

Philip Trettner's avatar
Philip Trettner committed
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
// intersection between point and obj is same as contains
template <int D, class ScalarT, class Obj, class = void_t<decltype(contains(std::declval<pos<D, ScalarT>>(), std::declval<Obj>()))>>
constexpr optional<pos<D, ScalarT>> intersection(pos<D, ScalarT> const& p, Obj const& obj)
{
    if (contains(obj, p))
        return p;
    return {};
}

// intersection between point and obj is same as contains
template <int D, class ScalarT, class Obj, class = void_t<decltype(contains(std::declval<pos<D, ScalarT>>(), std::declval<Obj>()))>>
constexpr optional<pos<D, ScalarT>> intersection(Obj const& obj, pos<D, ScalarT> const& p)
{
    if (contains(obj, p))
        return p;
    return {};
}


369
// ====================================== Ray Intersections from Line Intersections ======================================
370

371
372
template <int D, class ScalarT, class Obj>
[[nodiscard]] constexpr auto intersection_parameter(ray<D, ScalarT> const& ray, Obj const& obj)
373
    -> decltype(intersection_parameter(inf_of(ray), obj).is_hits, intersection_parameter(inf_of(ray), obj))
374
{
375
376
377
378
    const auto inter = intersection_parameter(inf_of(ray), obj);
    constexpr auto maxHits = inter.max_hits;
    if (!inter.any() || inter.last() < ScalarT(0))
        return {};
379

380
381
382
383
384
385
386
    if constexpr (maxHits == 2)
    {
        if (inter.first() < ScalarT(0))
            return {inter[1]};
    }
    else
        static_assert(maxHits == 1, "Only up to two intersections implemented");
387

388
    return inter;
389
}
390
391
template <int D, class ScalarT, class Obj>
[[nodiscard]] constexpr auto intersection_parameter(ray<D, ScalarT> const& ray, Obj const& obj)
392
    -> decltype(intersection_parameter(inf_of(ray), obj).value().is_hit_interval, intersection_parameter(inf_of(ray), obj))
393
{
394
395
396
397
398
399
400
401
402
403
404
405
406
    const auto inter = intersection_parameter(inf_of(ray), obj);

    if (!inter.has_value())
        return inter;

    auto interval = inter.value();
    if (interval.end < ScalarT(0))
        return {};

    TG_ASSERT((interval.start <= ScalarT(0)) == contains(obj, ray.origin));

    interval.start = max(interval.start, ScalarT(0));
    return interval;
407
408
}

409
// ====================================== Line - Object Intersections ======================================
Philip Trettner's avatar
Philip Trettner committed
410

411
// line - point
412
template <class ScalarT>
413
[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<1, ScalarT> const& l, pos<1, ScalarT> const& p)
414
{
415
    return coordinates(l, p);
416
417
}

418
// line - line
419
template <class ScalarT>
420
[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<2, ScalarT> const& l0, line<2, ScalarT> const& l1)
421
{
422
423
424
    // l0.pos + l0.dir * t.x == l1.pos + l1.dir * t.y  <=>  (l0.dir | -l1.dir) * (t.x | t.y)^T == l1.pos - l0.pos
    auto M = mat<2, 2, ScalarT>::from_cols(l0.dir, -l1.dir);
    auto t = inverse(M) * (l1.pos - l0.pos);
425
426
    if (!tg::is_finite(t.x))
        return {};
427
428
429
    return t.x;
}

430
// line - ray
Aaron Grabowy's avatar
Aaron Grabowy committed
431
template <class ScalarT>
432
[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<2, ScalarT> const& l, ray<2, ScalarT> const& r)
Aaron Grabowy's avatar
Aaron Grabowy committed
433
{
434
435
    auto M = mat<2, 2, ScalarT>::from_cols(l.dir, -r.dir);
    auto t = inverse(M) * (r.origin - l.pos);
436
    if (t.y < ScalarT(0) || !tg::is_finite(t.x))
Aaron Grabowy's avatar
Aaron Grabowy committed
437
438
439
440
        return {};
    return t.x;
}

441
// line - segment
442
template <class ScalarT>
443
[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<2, ScalarT> const& l, segment<2, ScalarT> const& s)
444
{
445
446
    auto M = mat<2, 2, ScalarT>::from_cols(l.dir, s.pos0 - s.pos1);
    auto t = inverse(M) * (s.pos0 - l.pos);
447
    if (t.y < ScalarT(0) || t.y > ScalarT(1) || !tg::is_finite(t.x))
448
449
450
451
        return {};
    return t.x;
}

452
// line - plane
453
template <int D, class ScalarT>
454
[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<D, ScalarT> const& l, plane<D, ScalarT> const& p)
455
{
456
457
    const auto dotND = dot(p.normal, l.dir);
    if (dotND == ScalarT(0)) // if plane normal and line direction are orthogonal, there is no intersection
Philip Trettner's avatar
Philip Trettner committed
458
        return {};
459

460
461
    // <l.pos + t * l.dir, p.normal> = p.dis  <=>  t = (p.dis - <l.pos, p.normal>) / <l.dir, p.normal>
    return (p.dis - dot(p.normal, l.pos)) / dotND;
462
463
}

464
// line - halfspace
465
template <int D, class ScalarT>
466
[[nodiscard]] constexpr optional<hit_interval<ScalarT>> intersection_parameter(line<D, ScalarT> const& l, halfspace<D, ScalarT> const& h)
467
{
468
469
470
471
    const auto dotND = dot(h.normal, l.dir);
    const auto dist = signed_distance(l.pos, h);

    if (dotND == ScalarT(0)) // if halfspace normal and line direction are orthogonal, there is no intersection
472
    {
473
474
475
        if (dist <= ScalarT(0))
            return {{tg::min<ScalarT>(), tg::max<ScalarT>()}}; // completely contained
        return {}; // completely outside
476
477
    }

478
479
480
481
    const auto t = -dist / dotND;
    if (dotND < ScalarT(0))
        return {{t, tg::max<ScalarT>()}}; // line goes into the halfspace
    return {{tg::min<ScalarT>(), t}}; // line goes out of the halfspace
482
483
484
}
template <int D, class ScalarT>
[[nodiscard]] constexpr optional<ScalarT> closest_intersection_parameter(ray<D, ScalarT> const& r, halfspace<D, ScalarT> const& h)
485
486
487
488
489
490
491
{
    // check if ray origin is already contained in the halfspace
    const auto dist = signed_distance(r.origin, h);
    if (dist <= ScalarT(0))
        return ScalarT(0);

    // if ray points away from the halfspace there is no intersection
492
    const auto dotND = dot(h.normal, r.dir);
493
494
495
    if (dotND >= ScalarT(0))
        return {};

496
    return -dist / dotND;
497
498
}

499
// line - aabb
Philip Trettner's avatar
Philip Trettner committed
500
template <int D, class ScalarT>
501
[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, aabb_boundary<D, ScalarT> const& b)
502
{
Aaron Grabowy's avatar
Aaron Grabowy committed
503
    // based on ideas from https://gamedev.stackexchange.com/q/18436
504
505
506
    auto tFirst = tg::min<ScalarT>();
    auto tSecond = tg::max<ScalarT>();
    for (auto i = 0; i < D; ++i)
507
    {
508
        if (abs(l.dir[i]) > ScalarT(100) * tg::epsilon<ScalarT>)
509
        {
510
511
            const auto tMin = (b.min[i] - l.pos[i]) / l.dir[i];
            const auto tMax = (b.max[i] - l.pos[i]) / l.dir[i];
512
513
514
            tFirst = max(tFirst, min(tMin, tMax));
            tSecond = min(tSecond, max(tMin, tMax));
        }
515
        else if (l.pos[i] < b.min[i] || l.pos[i] > b.max[i])
516
            return {}; // ray parallel to this axis and outside of the aabb
517
    }
Philip Trettner's avatar
Philip Trettner committed
518

519
520
    // no intersection if line misses the aabb
    if (tFirst > tSecond)
521
522
523
524
525
        return {};

    return {tFirst, tSecond};
}

526
// line - box
527
template <int D, class ScalarT>
528
[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, box_boundary<D, ScalarT> const& b)
529
{
530
531
    const auto bMin = b[comp<D, ScalarT>(-1)] - l.pos;
    const auto bMax = b[comp<D, ScalarT>(1)] - l.pos;
532
533
534
535
    auto tFirst = tg::min<ScalarT>();
    auto tSecond = tg::max<ScalarT>();
    for (auto i = 0; i < D; ++i)
    {
536
        const auto rDir = dot(l.dir, b.half_extents[i]);
537
538
539
540
541
542
543
544
545
546
547
        if (abs(rDir) > ScalarT(100) * tg::epsilon<ScalarT>)
        {
            const auto tMin = dot(bMin, b.half_extents[i]) / rDir;
            const auto tMax = dot(bMax, b.half_extents[i]) / rDir;
            tFirst = max(tFirst, min(tMin, tMax));
            tSecond = min(tSecond, max(tMin, tMax));
        }
        else if (dot(bMin, b.half_extents[i]) > ScalarT(0) || dot(bMax, b.half_extents[i]) < ScalarT(0))
            return {}; // ray parallel to this half_extents axis and outside of the box
    }

548
549
    // no intersection if line misses the box
    if (tFirst > tSecond)
550
        return {};
Philip Trettner's avatar
Philip Trettner committed
551

552
    return {tFirst, tSecond};
553
}
554
template <class ScalarT>
555
[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<3, ScalarT> const& l, box<2, ScalarT, 3> const& b)
556
{
557
558
    const auto t = intersection_parameter(l, plane<3, ScalarT>(normal_of(b), b.center));
    if (!t.any()) // line parallel to box plane
559
560
        return {};

561
    const auto p = l[t.first()] - b.center;
562
563
564
565
566
    if (abs(dot(b.half_extents[0], p)) > length_sqr(b.half_extents[0]) || abs(dot(b.half_extents[1], p)) > length_sqr(b.half_extents[1]))
        return {};

    return t;
}
567

568
// line - disk
569
template <class ScalarT>
570
[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<3, ScalarT> const& l, sphere<2, ScalarT, 3> const& d)
571
{
572
573
    const auto t = intersection_parameter(l, plane<3, ScalarT>(d.normal, d.center));
    if (!t.any()) // line parallel to disk plane
574
575
        return {};

576
    const auto p = l[t.first()];
577
578
579
580
581
582
    if (distance_sqr(p, d.center) > d.radius * d.radius)
        return {};

    return t;
}

583
// line - sphere
Philip Trettner's avatar
Philip Trettner committed
584
template <int D, class ScalarT>
585
[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, sphere_boundary<D, ScalarT> const& s)
586
{
587
    auto t = dot(s.center - l.pos, l.dir);
588

589
    auto d_sqr = distance_sqr(l[t], s.center);
Philip Trettner's avatar
Philip Trettner committed
590
591
592
593
594
    auto r_sqr = s.radius * s.radius;
    if (d_sqr > r_sqr)
        return {};

    auto dt = sqrt(r_sqr - d_sqr);
595
    return {t - dt, t + dt};
Philip Trettner's avatar
Philip Trettner committed
596
}
Philip Trettner's avatar
Philip Trettner committed
597

598
// line - hemisphere
599
template <int D, class ScalarT>
600
[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, hemisphere_boundary<D, ScalarT> const& h)
601
{
602
    return detail::merge_hits(l, caps_of(h), boundary_no_caps_of(h));
603
604
}
template <int D, class ScalarT>
605
[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, hemisphere_boundary_no_caps<D, ScalarT> const& h)
606
{
Aaron Grabowy's avatar
Aaron Grabowy committed
607
    ScalarT hits[2] = {};
608
    auto numHits = 0;
609
    const auto sphereHits = intersection_parameter(l, sphere_boundary<D, ScalarT>(h.center, h.radius));
610
611
    const auto halfSpace = halfspace<D, ScalarT>(-h.normal, h.center); // the intersection of this halfspace and the sphere is exactly the hemisphere
    for (const auto& hit : sphereHits)
612
        if (contains(halfSpace, l[hit]))
613
614
615
616
            hits[numHits++] = hit;

    return {hits, numHits};
}
617
618
template <class Obj, int D, class ScalarT, class TraitsT>
[[nodiscard]] constexpr bool intersects(Obj const& lr, hemisphere<D, ScalarT, TraitsT> const& h)
619
{
620
    static_assert(object_traits<Obj>::is_infinite, "For finite objects, complete containment within boundary has to be considered as well");
621
    if constexpr (std::is_same_v<TraitsT, boundary_no_caps_tag>)
622
        return intersection_parameter(lr, h).any();
623
    else
624
        return detail::intersects_any(lr, caps_of(h), boundary_no_caps_of(h));
625
626
}

627
// line - capsule
628
template <class ScalarT>
629
[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, capsule_boundary<3, ScalarT> const& c)
630
631
632
{
    using caps_t = hemisphere_boundary_no_caps<3, ScalarT>;
    const auto n = direction(c);
633
    return detail::merge_hits(l, caps_t(c.axis.pos0, c.radius, -n), caps_t(c.axis.pos1, c.radius, n),
634
635
                              cylinder_boundary_no_caps<3, ScalarT>(c.axis, c.radius));
}
636
637
template <class Obj, class ScalarT, class TraitsT>
[[nodiscard]] constexpr bool intersects(Obj const& lr, capsule<3, ScalarT, TraitsT> const& c)
638
{
639
    static_assert(object_traits<Obj>::is_infinite, "For finite objects, complete containment within boundary has to be considered as well");
640
    using caps_t = sphere_boundary<3, ScalarT>; // spheres are faster than hemispheres and equivalent for the yes/no decision
641
    return detail::intersects_any(lr, caps_t(c.axis.pos0, c.radius), caps_t(c.axis.pos1, c.radius),
642
643
644
                                  cylinder_boundary_no_caps<3, ScalarT>(c.axis, c.radius));
}

645
// line - cylinder
Philip Trettner's avatar
Philip Trettner committed
646
template <class ScalarT>
647
[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, cylinder_boundary<3, ScalarT> const& c)
Philip Trettner's avatar
Philip Trettner committed
648
{
649
    const auto caps = caps_of(c);
650
    return detail::merge_hits(l, caps[0], caps[1], boundary_no_caps_of(c));
651
}
652
653
template <class Obj, class ScalarT, class TraitsT>
[[nodiscard]] constexpr bool intersects(Obj const& lr, cylinder<3, ScalarT, TraitsT> const& c)
654
{
655
    static_assert(object_traits<Obj>::is_infinite, "For finite objects, complete containment within boundary has to be considered as well");
656
    if constexpr (std::is_same_v<TraitsT, boundary_no_caps_tag>)
657
        return intersection_parameter(lr, c).any();
658
659
660
    else
    {
        const auto caps = caps_of(c);
661
        return detail::intersects_any(lr, caps[0], caps[1], boundary_no_caps_of(c));
662
    }
663
}
664
template <class ScalarT>
665
[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, cylinder_boundary_no_caps<3, ScalarT> const& c)
666
{
667
668
669
    auto const infInter = intersection_parameter(l, inf_of(c));
    if (!infInter.any())
        return infInter;
670

671
672
673
    auto const d = c.axis.pos1 - c.axis.pos0;
    auto const lambda0 = dot(l[infInter[0]] - c.axis.pos0, d);
    auto const lambda1 = dot(l[infInter[1]] - c.axis.pos0, d);
674

Aaron Grabowy's avatar
Aaron Grabowy committed
675
    ScalarT hits[2] = {};
676
677
678
679
680
681
    auto numHits = 0;
    auto const dDotD = dot(d, d);
    if (ScalarT(0) <= lambda0 && lambda0 <= dDotD)
        hits[numHits++] = infInter[0];
    if (ScalarT(0) <= lambda1 && lambda1 <= dDotD)
        hits[numHits++] = infInter[1];
682

683
684
    return {hits, numHits};
}
685

686
687
688
689
690
691
// line - inf_cylinder
template <class ScalarT>
[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, inf_cylinder_boundary<3, ScalarT> const& c)
{
    auto const cosA = dot(c.axis.dir, l.dir);
    auto const sinA_sqr = 1 - cosA * cosA;
692

693
694
    if (sinA_sqr <= 0)
        return {}; // line and cylinder are parallel
695

696
697
698
699
700
701
    // compute closest points of line and cylinder axis
    auto const origDiff = l.pos - c.axis.pos;
    auto const fLine = dot(l.dir, origDiff);
    auto const fAxis = dot(c.axis.dir, origDiff);
    auto const tLine = (cosA * fAxis - fLine) / sinA_sqr;
    auto const tAxis = (fAxis - cosA * fLine) / sinA_sqr;
702

703
704
    auto const line_axis_dist_sqr = distance_sqr(l[tLine], c.axis[tAxis]);
    auto const cyl_radius_sqr = c.radius * c.radius;
705

706
707
    if (cyl_radius_sqr < line_axis_dist_sqr)
        return {}; // line misses the cylinder
708

709
710
711
712
713
    // radius in 2D slice = sqrt(cyl_radius_sqr - line_axis_dist_sqr)
    // infinite tube intersection
    auto const s = sqrt((cyl_radius_sqr - line_axis_dist_sqr) / sinA_sqr);
    return {tLine - s, tLine + s};
}
714
template <class ScalarT>
715
[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<2, ScalarT> const& l, inf_cylinder_boundary<2, ScalarT> const& c)
716
717
{
    const auto n = perpendicular(c.axis.dir);
718
719
    const auto d = dot(l.dir, n);
    if (d == ScalarT(0)) // line parallel to inf_cylinder
720
721
        return {};

722
    const auto dist = dot(c.axis.pos - l.pos, n);
723
    const auto [tMin, tMax] = minmax((dist - c.radius) / d, (dist + c.radius) / d);
724
    return {tMin, tMax};
725
726
}

Aaron Grabowy's avatar
Aaron Grabowy committed
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
// line - inf_cone
template <class ScalarT>
[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<2, ScalarT> const& l, inf_cone_boundary<2, ScalarT> const& c)
{
    auto ray1 = ray<2, ScalarT>(c.apex, rotate(c.opening_dir, c.opening_angle / ScalarT(2)));
    auto ray2 = ray<2, ScalarT>(c.apex, rotate(c.opening_dir, -c.opening_angle / ScalarT(2)));
    return detail::merge_hits(l, ray1, ray2);
}
template <class ScalarT>
[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, inf_cone_boundary<3, ScalarT> const& ic)
{
    // see https://lousodrome.net/blog/light/2017/01/03/intersection-of-a-ray-and-a-cone/
    auto const dv = dot(l.dir, ic.opening_dir);
    auto const cos2 = pow2(cos(ic.opening_angle * ScalarT(0.5)));
    auto const co = l.pos - ic.apex;
    auto const cov = dot(co, ic.opening_dir);
    auto const a = dv * dv - cos2;
    auto const b = ScalarT(2) * (dv * cov - dot(l.dir, co) * cos2);
    auto const c = cov * cov - dot(co, co) * cos2;
    auto const inter = detail::solve_quadratic(a, b, c);
Aaron Grabowy's avatar
Aaron Grabowy committed
747
748
    if (!inter.any())
        return inter;
Aaron Grabowy's avatar
Aaron Grabowy committed
749
750

    // exclude intersections with mirrored cone
Aaron Grabowy's avatar
Aaron Grabowy committed
751
    ScalarT hits[2] = {};
Aaron Grabowy's avatar
Aaron Grabowy committed
752
    auto numHits = 0;
Aaron Grabowy's avatar
Aaron Grabowy committed
753
    TG_ASSERT(ic.opening_angle <= tg::angle::from_degree(ScalarT(180)) && "Only convex objects are supported, but an inf_cone with openinge angle > 180 degree is not convex.");
Aaron Grabowy's avatar
Aaron Grabowy committed
754
755
756
757
    // if it is not used for solid cones, this works:
    // auto const coneDir = ic.opening_angle > 180_deg ? -ic.opening_dir : ic.opening_dir;
    // if (dot(l[inter[0]] - ic.apex, coneDir) >= ScalarT(0)) ...
    if (dot(l[inter[0]] - ic.apex, ic.opening_dir) >= ScalarT(0))
Aaron Grabowy's avatar
Aaron Grabowy committed
758
        hits[numHits++] = inter[0];
Aaron Grabowy's avatar
Aaron Grabowy committed
759
    if (dot(l[inter[1]] - ic.apex, ic.opening_dir) >= ScalarT(0))
Aaron Grabowy's avatar
Aaron Grabowy committed
760
761
762
763
764
        hits[numHits++] = inter[1];

    return {hits, numHits};
}

765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
// line - cone
template <class ScalarT>
[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, cone_boundary_no_caps<3, ScalarT> const& cone)
{
    auto const apex = apex_of(cone);
    auto const openingDir = -normal_of(cone.base);
    auto const borderPos = any_point(boundary_of(cone.base));
    auto const openingAngleHalf = angle_between(openingDir, normalize(borderPos - apex));

    // see https://lousodrome.net/blog/light/2017/01/03/intersection-of-a-ray-and-a-cone/
    auto const dv = dot(l.dir, openingDir);
    auto const cos2 = pow2(cos(openingAngleHalf));
    auto const co = l.pos - apex;
    auto const cov = dot(co, openingDir);
    auto const a = dv * dv - cos2;
    auto const b = ScalarT(2) * (dv * cov - dot(l.dir, co) * cos2);
    auto const c = cov * cov - dot(co, co) * cos2;
    auto const inter = detail::solve_quadratic(a, b, c);
    if (!inter.any())
        return inter;

    // exclude intersections with mirrored cone
Aaron Grabowy's avatar
Aaron Grabowy committed
787
    ScalarT hits[2] = {};
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
    auto numHits = 0;
    auto const h0 = dot(l[inter[0]] - apex, openingDir);
    auto const h1 = dot(l[inter[1]] - apex, openingDir);
    if (ScalarT(0) <= h0 && h0 <= cone.height)
        hits[numHits++] = inter[0];
    if (ScalarT(0) <= h1 && h1 <= cone.height)
        hits[numHits++] = inter[1];

    return {hits, numHits};
}

// line - pyramid
template <class BaseT, typename = std::enable_if_t<!std::is_same_v<BaseT, sphere<2, typename BaseT::scalar_t, 3>>>>
[[nodiscard]] constexpr hits<2, typename BaseT::scalar_t> intersection_parameter(line<3, typename BaseT::scalar_t> const& l, pyramid_boundary_no_caps<BaseT> const& py)
{
    auto const faces = faces_of(py);
Aaron Grabowy's avatar
Aaron Grabowy committed
804
    return detail::merge_hits_array(l, faces, std::make_index_sequence<faces.size()>{});
805
806
807
808
809
810
811
}
template <class BaseT>
[[nodiscard]] constexpr hits<2, typename BaseT::scalar_t> intersection_parameter(line<3, typename BaseT::scalar_t> const& l, pyramid_boundary<BaseT> const& py)
{
    return detail::merge_hits(l, py.base, boundary_no_caps_of(py));
}

812
// line - triangle2
813
template <class ScalarT>
814
[[nodiscard]] constexpr optional<hit_interval<ScalarT>> intersection_parameter(line<2, ScalarT> const& l, triangle<2, ScalarT> const& t)
815
816
{
    ScalarT closestIntersection = tg::max<ScalarT>();
817
    ScalarT furtherIntersection = tg::min<ScalarT>();
818
819
820
    auto numIntersections = 0;
    for (const auto& edge : edges_of(t))
    {