From 46ae55e7400e1fbe4b3231be743b76926090f1f6 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Wed, 4 Nov 2020 19:57:52 +0100
Subject: [PATCH 01/42] Added project(pos, pos)

---
 src/typed-geometry/functions/objects/project.hh | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/src/typed-geometry/functions/objects/project.hh b/src/typed-geometry/functions/objects/project.hh
index 5cd838e..a2f7857 100644
--- a/src/typed-geometry/functions/objects/project.hh
+++ b/src/typed-geometry/functions/objects/project.hh
@@ -60,6 +60,15 @@ template <int D, class ScalarT>
 }
 
 
+// ============== project to point ==============
+
+template <int D, class ScalarT>
+[[nodiscard]] constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, pos<D, ScalarT> const& q)
+{
+    return q;
+}
+
+
 // ============== project to line / ray / segment ==============
 
 template <int D, class ScalarT>
-- 
GitLab


From 693adea14919b8e60828ed92f8d0bc8fe6206735 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Wed, 4 Nov 2020 20:14:17 +0100
Subject: [PATCH 02/42] Fixed unused warning

---
 src/typed-geometry/functions/objects/project.hh | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/typed-geometry/functions/objects/project.hh b/src/typed-geometry/functions/objects/project.hh
index a2f7857..26747e8 100644
--- a/src/typed-geometry/functions/objects/project.hh
+++ b/src/typed-geometry/functions/objects/project.hh
@@ -63,7 +63,7 @@ template <int D, class ScalarT>
 // ============== project to point ==============
 
 template <int D, class ScalarT>
-[[nodiscard]] constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, pos<D, ScalarT> const& q)
+[[nodiscard]] constexpr pos<D, ScalarT> project([[maybe_unused]] pos<D, ScalarT> const& p, pos<D, ScalarT> const& q)
 {
     return q;
 }
-- 
GitLab


From 1f9e750a97130832148258738e93d706cc341c4f Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Fri, 6 Nov 2020 00:51:29 +0100
Subject: [PATCH 03/42] Reordered and grouped intersection functions

---
 .../functions/objects/intersection.hh         | 225 +++++++++---------
 1 file changed, 116 insertions(+), 109 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 98425f5..5904783 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -385,6 +385,111 @@ template <class ScalarT>
     return t;
 }
 
+// ray - triangle
+template <class ScalarT>
+[[nodiscard]] constexpr optional<ScalarT> intersection_parameter(ray<3, ScalarT> const& r,
+                                                                 triangle<3, ScalarT> const& t,
+                                                                 dont_deduce<ScalarT> eps = 100 * tg::epsilon<ScalarT>)
+{
+    auto e1 = t.pos1 - t.pos0;
+    auto e2 = t.pos2 - t.pos0;
+
+    auto pvec = tg::cross(r.dir, e2);
+    auto det = dot(pvec, e1);
+
+    if (det < ScalarT(0))
+    {
+        std::swap(e1, e2);
+        pvec = tg::cross(r.dir, e2);
+        det = -det;
+    }
+
+    if (det < eps)
+        return {};
+
+    auto tvec = r.origin - t.pos0;
+    auto u = dot(tvec, pvec);
+    if (u < ScalarT(0) || u > det)
+        return {};
+
+    auto qvec = cross(tvec, e1);
+    auto v = dot(r.dir, qvec);
+    if (v < ScalarT(0) || v + u > det)
+        return {};
+
+    auto lambda = (ScalarT(1) / det) * dot(e2, qvec);
+    return (lambda > ScalarT(0)) ? lambda : tg::optional<ScalarT>();
+}
+
+// ray - box
+template <class ScalarT>
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, box<3, ScalarT> const& b)
+{
+    // see https://github.com/gszauer/GamePhysicsCookbook/blob/master/Code/Geometry3D.cpp
+
+    auto const p = b.center - r.origin;
+
+    // TODO: without normalize
+    auto const X = normalize(b.half_extents[0]);
+    auto const Y = normalize(b.half_extents[1]);
+    auto const Z = normalize(b.half_extents[2]);
+
+    auto f = tg::vec3(dot(X, r.dir), //
+                      dot(Y, r.dir), //
+                      dot(Z, r.dir));
+
+    auto const e = tg::vec3(dot(X, p), //
+                            dot(Y, p), //
+                            dot(Z, p));
+
+    auto const size = tg::vec3(length(b.half_extents[0]), //
+                               length(b.half_extents[1]), //
+                               length(b.half_extents[2]));
+
+    float t[6] = {0, 0, 0, 0, 0, 0};
+    for (int i = 0; i < 3; ++i)
+    {
+        if (abs(f[i]) < tg::epsilon<ScalarT> / 1000)
+        {
+            if (-e[i] - size[i] > 0 || -e[i] + size[i] < 0)
+                return {};
+
+            f[i] = tg::epsilon<ScalarT>; // Avoid div by 0!
+        }
+
+        t[i * 2 + 0] = (e[i] + size[i]) / f[i]; // tmin[x, y, z]
+        t[i * 2 + 1] = (e[i] - size[i]) / f[i]; // tmax[x, y, z]
+    }
+
+    auto const tmin = tg::max(tg::max(tg::min(t[0], t[1]), tg::min(t[2], t[3])), tg::min(t[4], t[5]));
+    auto const tmax = tg::min(tg::min(tg::max(t[0], t[1]), tg::max(t[2], t[3])), tg::max(t[4], t[5]));
+
+    // if tmax < 0, ray is intersecting AABB
+    // but entire AABB is behing it's origin
+    if (tmax < 0)
+        return {};
+
+    // if tmin > tmax, ray doesn't intersect AABB
+    if (tmin > tmax)
+        return {};
+
+    ScalarT hits[2];
+
+    if (tmin >= 0) // two valid intersections
+    {
+        hits[0] = tmin;
+        hits[1] = tmax;
+        return {hits, 2};
+    }
+    else // one valid intersection
+    {
+        hits[0] = tmax;
+        return {hits, 1};
+    }
+}
+
+
+// ====================================== Object - Object Intersections ======================================
 
 // returns intersection circle of sphere and sphere (normal points from a to b)
 // for now does not work if spheres are identical (result would be a sphere3 again)
@@ -613,106 +718,23 @@ template <int D, class ScalarT>
     return t;
 }
 
-template <class ScalarT>
-[[nodiscard]] constexpr optional<ScalarT> intersection_parameter(ray<3, ScalarT> const& r,
-                                                                 triangle<3, ScalarT> const& t,
-                                                                 dont_deduce<ScalarT> eps = 100 * tg::epsilon<ScalarT>)
+template <int D, class ScalarT>
+[[nodiscard]] constexpr optional<segment<D, ScalarT>> intersection(segment<D, ScalarT> const& a, sphere<D, ScalarT> const& b)
 {
-    auto e1 = t.pos1 - t.pos0;
-    auto e2 = t.pos2 - t.pos0;
-
-    auto pvec = tg::cross(r.dir, e2);
-    auto det = dot(pvec, e1);
-
-    if (det < ScalarT(0))
-    {
-        std::swap(e1, e2);
-        pvec = tg::cross(r.dir, e2);
-        det = -det;
-    }
-
-    if (det < eps)
-        return {};
-
-    auto tvec = r.origin - t.pos0;
-    auto u = dot(tvec, pvec);
-    if (u < ScalarT(0) || u > det)
-        return {};
-
-    auto qvec = cross(tvec, e1);
-    auto v = dot(r.dir, qvec);
-    if (v < ScalarT(0) || v + u > det)
-        return {};
-
-    auto lambda = (ScalarT(1) / det) * dot(e2, qvec);
-    return (lambda > ScalarT(0)) ? lambda : tg::optional<ScalarT>();
+    static_assert(always_false<ScalarT>, "not implemented");
+    (void)a;
+    (void)b;
+    return {}; // TODO
 }
 
-template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, box<3, ScalarT> const& b)
+template <int D, class ScalarT>
+[[nodiscard]] constexpr optional<segment<D, ScalarT>> intersection(sphere<D, ScalarT> const& b, segment<D, ScalarT> const& a)
 {
-    // see https://github.com/gszauer/GamePhysicsCookbook/blob/master/Code/Geometry3D.cpp
-
-    auto const p = b.center - r.origin;
-
-    // TODO: without normalize
-    auto const X = normalize(b.half_extents[0]);
-    auto const Y = normalize(b.half_extents[1]);
-    auto const Z = normalize(b.half_extents[2]);
-
-    auto f = tg::vec3(dot(X, r.dir), //
-                      dot(Y, r.dir), //
-                      dot(Z, r.dir));
-
-    auto const e = tg::vec3(dot(X, p), //
-                            dot(Y, p), //
-                            dot(Z, p));
-
-    auto const size = tg::vec3(length(b.half_extents[0]), //
-                               length(b.half_extents[1]), //
-                               length(b.half_extents[2]));
-
-    float t[6] = {0, 0, 0, 0, 0, 0};
-    for (int i = 0; i < 3; ++i)
-    {
-        if (abs(f[i]) < tg::epsilon<ScalarT> / 1000)
-        {
-            if (-e[i] - size[i] > 0 || -e[i] + size[i] < 0)
-                return {};
-
-            f[i] = tg::epsilon<ScalarT>; // Avoid div by 0!
-        }
-
-        t[i * 2 + 0] = (e[i] + size[i]) / f[i]; // tmin[x, y, z]
-        t[i * 2 + 1] = (e[i] - size[i]) / f[i]; // tmax[x, y, z]
-    }
-
-    auto const tmin = tg::max(tg::max(tg::min(t[0], t[1]), tg::min(t[2], t[3])), tg::min(t[4], t[5]));
-    auto const tmax = tg::min(tg::min(tg::max(t[0], t[1]), tg::max(t[2], t[3])), tg::max(t[4], t[5]));
-
-    // if tmax < 0, ray is intersecting AABB
-    // but entire AABB is behing it's origin
-    if (tmax < 0)
-        return {};
+    return intersection(b, a);
+}
 
-    // if tmin > tmax, ray doesn't intersect AABB
-    if (tmin > tmax)
-        return {};
 
-    ScalarT hits[2];
-
-    if (tmin >= 0) // two valid intersections
-    {
-        hits[0] = tmin;
-        hits[1] = tmax;
-        return {hits, 2};
-    }
-    else // one valid intersection
-    {
-        hits[0] = tmax;
-        return {hits, 1};
-    }
-}
+// ====================================== Checks if Object Intersects aabb ======================================
 
 template <int D, class ScalarT>
 [[nodiscard]] constexpr bool intersects(sphere<D, ScalarT> const& a, aabb<D, ScalarT> const& b)
@@ -759,21 +781,6 @@ template <int D, class ScalarT>
     return intersects(b, a);
 }
 
-template <int D, class ScalarT>
-[[nodiscard]] constexpr optional<segment<D, ScalarT>> intersection(segment<D, ScalarT> const& a, sphere<D, ScalarT> const& b)
-{
-    static_assert(always_false<ScalarT>, "not implemented");
-    (void)a;
-    (void)b;
-    return {}; // TODO
-}
-
-template <int D, class ScalarT>
-[[nodiscard]] constexpr optional<segment<D, ScalarT>> intersection(sphere<D, ScalarT> const& b, segment<D, ScalarT> const& a)
-{
-    return intersection(b, a);
-}
-
 template <class ScalarT>
 [[nodiscard]] constexpr bool intersects(triangle<2, ScalarT> const& a, aabb<2, ScalarT> const& b)
 {
-- 
GitLab


From 089bef8dadb82bb3f1bfd424b03fcb47d7882959 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Fri, 6 Nov 2020 01:22:02 +0100
Subject: [PATCH 04/42] Some slight code simplifications

---
 .../functions/objects/intersection.hh         | 37 +++++++------------
 1 file changed, 13 insertions(+), 24 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 5904783..a19a5f9 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -191,7 +191,7 @@ template <int D, class ScalarT>
     // <r.origin + t * r.dir, p.normal> = p.dis
     // t = (p.dis - <r.origin, p.normal>) / <r.dir, p.normal>
 
-    auto t = (p.dis - dot(p.normal, vec<D, ScalarT>(r.origin))) / dotND;
+    auto t = (p.dis - dot(p.normal, r.origin)) / dotND;
 
     // check whether plane lies behind ray
     if (t < 0)
@@ -364,10 +364,10 @@ template <class ScalarT>
 template <class ScalarT>
 [[nodiscard]] constexpr optional<ScalarT> closest_intersection_parameter(ray<3, ScalarT> const& r, cylinder<3, ScalarT> const& c)
 {
-    auto const dir = direction(c);
-    auto const t_cyl = closest_intersection_parameter(r, cylinder_boundary_no_caps<3, ScalarT>(c.axis, c.radius));
-    auto const t_cap0 = intersection_parameter(r, sphere<2, ScalarT, 3>(c.axis.pos0, c.radius, dir));
-    auto const t_cap1 = intersection_parameter(r, sphere<2, ScalarT, 3>(c.axis.pos1, c.radius, dir));
+    const auto caps = caps_of(c);
+    const auto t_cap0 = intersection_parameter(r, caps[0]);
+    const auto t_cap1 = intersection_parameter(r, caps[1]);
+    const auto t_cyl = closest_intersection_parameter(r, boundary_no_caps_of(c));
 
     optional<ScalarT> t;
 
@@ -696,7 +696,7 @@ template <int D, class ScalarT>
 template <class ScalarT>
 [[nodiscard]] constexpr ScalarT intersection_parameter(line<3, ScalarT> const& l, plane<3, ScalarT> const& p)
 {
-    return (p.dis - dot(l.pos - pos<3, ScalarT>::zero, p.normal)) / dot(l.dir, p.normal);
+    return (p.dis - dot(l.pos, p.normal)) / dot(l.dir, p.normal);
 }
 
 template <class ScalarT>
@@ -712,7 +712,7 @@ template <int D, class ScalarT>
     if (denom == ScalarT(0))
         return {};
 
-    auto t = (p.dis - dot(p.normal, a.pos0 - tg::pos<D, ScalarT>::zero)) / denom;
+    auto t = (p.dis - dot(p.normal, a.pos0)) / denom;
     if (t < 0 || t > 1)
         return {};
     return t;
@@ -784,23 +784,16 @@ template <int D, class ScalarT>
 template <class ScalarT>
 [[nodiscard]] constexpr bool intersects(triangle<2, ScalarT> const& a, aabb<2, ScalarT> const& b)
 {
+    if (!intersects(aabb_of(a), b))
+        return false;
+
     auto p0 = a.pos0;
     auto p1 = a.pos1;
     auto p2 = a.pos2;
-
-    auto bb_t = aabb_of(p0, p1, p2);
-    if (!intersects(bb_t, b))
-        return false;
-
     if (contains(b, p0) || contains(b, p1) || contains(b, p2))
         return true;
 
-    pos<2, ScalarT> aabb_pts[] = {
-        {b.min.x, b.min.y}, //
-        {b.min.x, b.max.y}, //
-        {b.max.x, b.min.y}, //
-        {b.max.x, b.max.y}, //
-    };
+    auto aabb_pts = vertices_of(b);
 
     auto const is_separate = [&](pos<2, ScalarT> pa, vec<2, ScalarT> n, pos<2, ScalarT> pb) {
         auto da = dot(n, pa);
@@ -847,7 +840,7 @@ template <class ScalarT>
     using pos_t = pos<3, ScalarT>;
     using vec_t = vec<3, ScalarT>;
 
-    auto const center = (bb_in.max + bb_in.min) / ScalarT(2);
+    auto const center = centroid_of(bb_in);
     auto const amin = pos_t(bb_in.min - center);
     auto const amax = pos_t(bb_in.max - center);
     auto const bb = aabb<3, ScalarT>(amin, amax);
@@ -868,12 +861,8 @@ template <class ScalarT>
                b.min.z < p.z && p.z < b.max.z;
     };
 
-    auto const contains_p0 = proper_contains(bb, p0);
-    auto const contains_p1 = proper_contains(bb, p1);
-    auto const contains_p2 = proper_contains(bb, p2);
-
     // early in: tri points vs AABB
-    if (contains_p0 || contains_p1 || contains_p2)
+    if (proper_contains(bb, p0) || proper_contains(bb, p1) || proper_contains(bb, p2))
         return true;
 
     // get adjusted tri base plane
-- 
GitLab


From 82ee1b49a081c6b62cb6f63a4e7d6d4edbcadc3e Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Fri, 6 Nov 2020 13:52:14 +0100
Subject: [PATCH 05/42] Moved cylinder_boundary_no_caps intersection next to
 cylinder intersection

---
 .../functions/objects/intersection.hh         | 100 +++++++++---------
 1 file changed, 50 insertions(+), 50 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index a19a5f9..6fa73ac 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -252,56 +252,6 @@ template <int D, class ScalarT>
     return {};
 }
 
-// ray - tube
-template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, cylinder_boundary_no_caps<3, ScalarT> const& c)
-{
-    auto cdir = direction(c);
-    auto cosA = dot(cdir, r.dir);
-    auto sinA_sqr = 1 - cosA * cosA;
-
-    if (sinA_sqr <= 0)
-        return {};
-
-    // compute closest points of the two lines
-    auto origDiff = r.origin - c.axis.pos0;
-    auto fRay = dot(r.dir, origDiff);
-    auto fLine = dot(cdir, origDiff);
-    auto tRay = (cosA * fLine - fRay) / sinA_sqr;
-    auto tLine = (fLine - cosA * fRay) / sinA_sqr;
-
-    auto closest_on_ray = r.origin + tRay * r.dir;
-    auto closest_on_line = c.axis.pos0 + tLine * cdir;
-    auto line_ray_dist_sqr = distance_sqr(closest_on_ray, closest_on_line);
-    auto cyl_radius_sqr = c.radius * c.radius;
-
-    if (line_ray_dist_sqr > cyl_radius_sqr)
-        return {};
-
-    // radius in 2D slice
-    auto r_2D = sqrt(cyl_radius_sqr - line_ray_dist_sqr);
-
-    // infinite tube intersection
-    auto s = r_2D / sqrt(sinA_sqr);
-    auto cyl_intersection_0 = closest_on_ray - s * r.dir;
-    auto cyl_intersection_1 = closest_on_ray + s * r.dir;
-
-    // project onto line segment
-    auto line_length = length(c.axis);
-    auto lambda_0 = dot(cyl_intersection_0 - c.axis.pos0, cdir);
-    auto lambda_1 = dot(cyl_intersection_1 - c.axis.pos0, cdir);
-
-    ScalarT hits[2];
-    int hit_cnt = 0;
-
-    if (tRay - s >= 0 && 0 <= lambda_0 && lambda_0 < line_length)
-        hits[hit_cnt++] = tRay - s;
-    if (tRay + s >= 0 && 0 <= lambda_1 && lambda_1 < line_length)
-        hits[hit_cnt++] = tRay + s;
-
-    return {hits, hit_cnt};
-}
-
 // ray - disk
 template <class ScalarT>
 [[nodiscard]] constexpr optional<ScalarT> intersection_parameter(ray<3, ScalarT> const& r, sphere<2, ScalarT, 3> const& d)
@@ -385,6 +335,56 @@ template <class ScalarT>
     return t;
 }
 
+// ray - tube
+template <class ScalarT>
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, cylinder_boundary_no_caps<3, ScalarT> const& c)
+{
+    auto cdir = direction(c);
+    auto cosA = dot(cdir, r.dir);
+    auto sinA_sqr = 1 - cosA * cosA;
+
+    if (sinA_sqr <= 0)
+        return {};
+
+    // compute closest points of the two lines
+    auto origDiff = r.origin - c.axis.pos0;
+    auto fRay = dot(r.dir, origDiff);
+    auto fLine = dot(cdir, origDiff);
+    auto tRay = (cosA * fLine - fRay) / sinA_sqr;
+    auto tLine = (fLine - cosA * fRay) / sinA_sqr;
+
+    auto closest_on_ray = r.origin + tRay * r.dir;
+    auto closest_on_line = c.axis.pos0 + tLine * cdir;
+    auto line_ray_dist_sqr = distance_sqr(closest_on_ray, closest_on_line);
+    auto cyl_radius_sqr = c.radius * c.radius;
+
+    if (line_ray_dist_sqr > cyl_radius_sqr)
+        return {};
+
+    // radius in 2D slice
+    auto r_2D = sqrt(cyl_radius_sqr - line_ray_dist_sqr);
+
+    // infinite tube intersection
+    auto s = r_2D / sqrt(sinA_sqr);
+    auto cyl_intersection_0 = closest_on_ray - s * r.dir;
+    auto cyl_intersection_1 = closest_on_ray + s * r.dir;
+
+    // project onto line segment
+    auto line_length = length(c.axis);
+    auto lambda_0 = dot(cyl_intersection_0 - c.axis.pos0, cdir);
+    auto lambda_1 = dot(cyl_intersection_1 - c.axis.pos0, cdir);
+
+    ScalarT hits[2];
+    int hit_cnt = 0;
+
+    if (tRay - s >= 0 && 0 <= lambda_0 && lambda_0 < line_length)
+        hits[hit_cnt++] = tRay - s;
+    if (tRay + s >= 0 && 0 <= lambda_1 && lambda_1 < line_length)
+        hits[hit_cnt++] = tRay + s;
+
+    return {hits, hit_cnt};
+}
+
 // ray - triangle
 template <class ScalarT>
 [[nodiscard]] constexpr optional<ScalarT> intersection_parameter(ray<3, ScalarT> const& r,
-- 
GitLab


From b8342f1998704bec1e13a4902228493560b19c88 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Sat, 7 Nov 2020 15:06:50 +0100
Subject: [PATCH 06/42] Added variadic constructor to ray_hits; Changed all
 intersection_parameter(ray, obj) to return ray_hits

---
 .../functions/objects/intersection.hh         | 104 +++++++++---------
 1 file changed, 49 insertions(+), 55 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 6fa73ac..aac22b2 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -27,16 +27,16 @@
 // intersects(a, b)              -> bool
 // intersection(a, b)            -> ???
 // intersection_safe(a, b)       -> optional<???>
-// intersection_parameter(a, b)  -> coords?
+// intersection_parameter(a, b)  -> coords? (for a ray: ray_hits<N, ScalarT>)
 // intersection_parameters(a, b) -> pair<coords, coords>?
 // intersection_exact(a, b)      -> variant
-// closest_intersection(a, b)            -> ???
-// closest_intersection_parameter(a, b)  -> coords?
+// closest_intersection(a, b)            -> position/object (for a ray: optional<pos>)
+// closest_intersection_parameter(a, b)  -> coords (for a ray: optional<ScalarT>)
 
 // "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
-// "intersection_parameter" returns a coordinates for the first object such that a[coords] == intersection(a, b)
+// "intersection_parameter" returns coordinates for the first object such that a[coords] == intersection(a, b)
 // "intersection_parameters" returns coordinates for both objects
 // "intersection_exact" returns a variant type describing all possible intersections, including degenerate cases
 // the "closest_" variants only return the closest intersection for objects where that concept is applicable (e.g. for rays)
@@ -46,6 +46,11 @@
 //  - intersection_safe is currently unsupported
 //  - for more elaborate ray-tracing, a future ray_cast function will exist (which also returns the intersection normal)
 
+// Implementation guidelines:
+// explicit intersection_parameter(ray, obj), which gives closest_intersection_parameter and intersection position
+// explicit closest_intersection_parameter(ray, obj) if this is faster than computing all intersections
+// explicit intersects(obj, aabb), which gives intersects(aabb, obj)
+
 
 namespace tg
 {
@@ -65,17 +70,22 @@ struct ray_hits
     template <class OtherT>
     using as_ray_hits = ray_hits<MaxHits, OtherT>;
 
-    int size() const { return _size; }
-    bool any() const { return _size > 0; }
+    [[nodiscard]] int size() const { return _size; }
+    [[nodiscard]] bool any() const { return _size > 0; }
 
-    HitT const& operator[](int idx)
+    HitT const& operator[](int idx) const
     {
         TG_ASSERT(0 <= idx && idx < _size);
         return _hit[idx];
     }
+    [[nodiscard]] HitT const& first() const
+    {
+        TG_ASSERT(_size > 0);
+        return _hit[0];
+    }
 
-    HitT const* begin() const { return _hit; }
-    HitT const* end() const { return _hit + _size; }
+    [[nodiscard]] HitT const* begin() const { return _hit; }
+    [[nodiscard]] HitT const* end() const { return _hit + _size; }
 
     ray_hits() = default;
     ray_hits(HitT* hits, int size) : _size(size)
@@ -83,6 +93,10 @@ struct ray_hits
         for (auto i = 0; i < size; ++i)
             _hit[i] = hits[i];
     }
+    template <typename... HitTs>
+    ray_hits(HitTs... hits) : _size(sizeof...(HitTs)), _hit{hits...}
+    {
+    }
 
 private:
     int _size = 0;
@@ -100,6 +114,13 @@ template <class A, class B>
     return intersection(a, b).has_value();
 }
 
+// if ray_hits intersection parameter is available and applicable, use that
+template <int D, class ScalarT, class Obj>
+[[nodiscard]] constexpr auto intersects(ray<D, ScalarT> const& r, Obj const& obj) -> decltype(intersection_parameter(r, obj).any())
+{
+    return intersection_parameter(r, obj).any();
+}
+
 // if a value-typed intersection parameter is available and applicable, use that
 template <class A, class B>
 [[nodiscard]] constexpr auto intersection(A const& a, B const& b) -> decltype(a[intersection_parameter(a, b)])
@@ -151,7 +172,7 @@ template <int D, class ScalarT, class Obj>
 {
     auto hits = intersection_parameter(r, obj);
     if (hits.any())
-        return hits[0];
+        return hits.first();
     return {};
 }
 
@@ -178,7 +199,7 @@ constexpr optional<pos<D, ScalarT>> intersection(Obj const& obj, pos<D, ScalarT>
 
 // ray - plane
 template <int D, class ScalarT>
-[[nodiscard]] constexpr optional<ScalarT> intersection_parameter(ray<D, ScalarT> const& r, plane<D, ScalarT> const& p)
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, plane<D, ScalarT> const& p)
 {
     // if plane normal and ray direction are parallel there is no intersection
     auto dotND = dot(p.normal, r.dir);
@@ -194,7 +215,7 @@ template <int D, class ScalarT>
     auto t = (p.dis - dot(p.normal, r.origin)) / dotND;
 
     // check whether plane lies behind ray
-    if (t < 0)
+    if (t < ScalarT(0))
         return {};
 
     return t;
@@ -234,33 +255,24 @@ template <int D, class ScalarT>
 
     auto dt = sqrt(r_sqr - d_sqr);
 
-    ScalarT hits[2];
-
     if (t - dt >= 0)
-    {
-        hits[0] = t - dt;
-        hits[1] = t + dt;
-        return {hits, 2};
-    }
+        return {t - dt, t + dt};
 
     if (t + dt >= 0)
-    {
-        hits[0] = t + dt;
-        return {hits, 1};
-    }
+        return {t + dt};
 
     return {};
 }
 
 // ray - disk
 template <class ScalarT>
-[[nodiscard]] constexpr optional<ScalarT> intersection_parameter(ray<3, ScalarT> const& r, sphere<2, ScalarT, 3> const& d)
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, sphere<2, ScalarT, 3> const& d)
 {
     auto const t = intersection_parameter(r, plane<3, ScalarT>(d.normal, d.center));
-    if (!t.has_value())
+    if (!t.any())
         return {};
 
-    auto const p = r[t.value()];
+    auto const p = r[t.first()];
     if (distance_sqr(p, d.center) > d.radius * d.radius)
         return {};
 
@@ -292,27 +304,18 @@ template <class ScalarT>
     auto tMin = min(t1, t2);
     auto tMax = max(t1, t2);
 
-    ScalarT hits[2];
-
     if (tMin >= ScalarT(0))
-    {
-        hits[0] = tMin;
-        hits[1] = tMax;
-        return {hits, 2};
-    }
+        return {tMin, tMax};
 
     if (tMax >= ScalarT(0))
-    {
-        hits[0] = tMax;
-        return {hits, 1};
-    }
+        return tMax;
 
     return {};
 }
 
 // ray - cylinder
 template <class ScalarT>
-[[nodiscard]] constexpr optional<ScalarT> closest_intersection_parameter(ray<3, ScalarT> const& r, cylinder<3, ScalarT> const& c)
+[[nodiscard]] constexpr optional<ScalarT> closest_intersection_parameter(ray<3, ScalarT> const& r, cylinder_boundary<3, ScalarT> const& c)
 {
     const auto caps = caps_of(c);
     const auto t_cap0 = intersection_parameter(r, caps[0]);
@@ -324,11 +327,11 @@ template <class ScalarT>
     if (t_cyl.has_value())
         t = t_cyl.value();
 
-    if (t_cap0.has_value() && (!t.has_value() || t.value() > t_cap0.value()))
-        t = t_cap0;
+    if (t_cap0.any() && (!t.has_value() || t.value() > t_cap0.first()))
+        t = t_cap0.first();
 
-    if (t_cap1.has_value() && (!t.has_value() || t.value() > t_cap1.value()))
-        t = t_cap1;
+    if (t_cap1.any() && (!t.has_value() || t.value() > t_cap1.first()))
+        t = t_cap1.first();
 
     TG_INTERNAL_ASSERT(!t.has_value() || t.value() >= 0);
 
@@ -387,7 +390,7 @@ template <class ScalarT>
 
 // ray - triangle
 template <class ScalarT>
-[[nodiscard]] constexpr optional<ScalarT> intersection_parameter(ray<3, ScalarT> const& r,
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<3, ScalarT> const& r,
                                                                  triangle<3, ScalarT> const& t,
                                                                  dont_deduce<ScalarT> eps = 100 * tg::epsilon<ScalarT>)
 {
@@ -418,7 +421,7 @@ template <class ScalarT>
         return {};
 
     auto lambda = (ScalarT(1) / det) * dot(e2, qvec);
-    return (lambda > ScalarT(0)) ? lambda : tg::optional<ScalarT>();
+    return (lambda > ScalarT(0)) ? lambda : ray_hits<1, ScalarT>();
 }
 
 // ray - box
@@ -473,19 +476,10 @@ template <class ScalarT>
     if (tmin > tmax)
         return {};
 
-    ScalarT hits[2];
-
     if (tmin >= 0) // two valid intersections
-    {
-        hits[0] = tmin;
-        hits[1] = tmax;
-        return {hits, 2};
-    }
+        return {tmin, tmax};
     else // one valid intersection
-    {
-        hits[0] = tmax;
-        return {hits, 1};
-    }
+        return tmax;
 }
 
 
-- 
GitLab


From 2d80567b1f6538d1f8e93e2e75dc7716a2984fe6 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Sat, 7 Nov 2020 15:16:37 +0100
Subject: [PATCH 07/42] Added ray2 - ray2 intersection

---
 src/typed-geometry/functions/objects/intersection.hh | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index aac22b2..82c3c73 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -197,6 +197,17 @@ constexpr optional<pos<D, ScalarT>> intersection(Obj const& obj, pos<D, ScalarT>
 
 // ====================================== Ray - Object Intersections ======================================
 
+// ray - ray
+template <class ScalarT>
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<2, ScalarT> const& r0, ray<2, ScalarT> const& r1)
+{
+    auto M = mat<2, 2, ScalarT>::from_cols(r0.dir, -r1.dir);
+    auto t = inverse(M) * (r1.origin - r0.origin);
+    if (t.x < ScalarT(0) || t.y < ScalarT(0))
+        return {};
+    return t.x;
+}
+
 // ray - plane
 template <int D, class ScalarT>
 [[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, plane<D, ScalarT> const& p)
-- 
GitLab


From 1c62c5d3f03d23aead25bdbc35d5e46f38cc8f6d Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Sat, 7 Nov 2020 15:48:21 +0100
Subject: [PATCH 08/42] Added ray2 - line2 and ray2 - segment2 intersections

---
 .../functions/objects/intersection.hh         | 22 +++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 82c3c73..1fd6855 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -197,6 +197,17 @@ constexpr optional<pos<D, ScalarT>> intersection(Obj const& obj, pos<D, ScalarT>
 
 // ====================================== Ray - Object Intersections ======================================
 
+// ray - line
+template <class ScalarT>
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<2, ScalarT> const& r, line<2, ScalarT> const& l)
+{
+    auto M = mat<2, 2, ScalarT>::from_cols(r.dir, -l.dir);
+    auto t = inverse(M) * (l.pos - r.origin);
+    if (t.x < ScalarT(0))
+        return {};
+    return t.x;
+}
+
 // ray - ray
 template <class ScalarT>
 [[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<2, ScalarT> const& r0, ray<2, ScalarT> const& r1)
@@ -208,6 +219,17 @@ template <class ScalarT>
     return t.x;
 }
 
+// ray - segment
+template <class ScalarT>
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<2, ScalarT> const& r, segment<2, ScalarT> const& s)
+{
+    auto M = mat<2, 2, ScalarT>::from_cols(r.dir, s.pos0 - s.pos1);
+    auto t = inverse(M) * (s.pos0 - r.origin);
+    if (t.x < ScalarT(0) || t.y < ScalarT(0) || t.y > ScalarT(1))
+        return {};
+    return t.x;
+}
+
 // ray - plane
 template <int D, class ScalarT>
 [[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, plane<D, ScalarT> const& p)
-- 
GitLab


From c93c3f4435187d8a012d5217f6342cc9205663de Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Sun, 8 Nov 2020 22:58:17 +0100
Subject: [PATCH 09/42] Default intersects(aabb, obj) implementation with
 reversed parameter order

---
 .../functions/objects/intersection.hh         | 30 +++++++++----------
 1 file changed, 14 insertions(+), 16 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 1fd6855..c8e0544 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -47,7 +47,7 @@
 //  - for more elaborate ray-tracing, a future ray_cast function will exist (which also returns the intersection normal)
 
 // Implementation guidelines:
-// explicit intersection_parameter(ray, obj), which gives closest_intersection_parameter and intersection position
+// explicit intersection_parameter(ray, obj), which gives closest_intersection_parameter and intersection position(s)
 // explicit closest_intersection_parameter(ray, obj) if this is faster than computing all intersections
 // explicit intersects(obj, aabb), which gives intersects(aabb, obj)
 
@@ -121,6 +121,19 @@ template <int D, class ScalarT, class Obj>
     return intersection_parameter(r, obj).any();
 }
 
+// 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();
+}
+
 // if a value-typed intersection parameter is available and applicable, use that
 template <class A, class B>
 [[nodiscard]] constexpr auto intersection(A const& a, B const& b) -> decltype(a[intersection_parameter(a, b)])
@@ -802,11 +815,6 @@ template <int D, class ScalarT>
 
     return d_min <= a.radius * a.radius;
 }
-template <int D, class ScalarT>
-[[nodiscard]] constexpr bool intersects(aabb<D, ScalarT> const& a, sphere<D, ScalarT> const& b)
-{
-    return intersects(b, a);
-}
 
 template <class ScalarT>
 [[nodiscard]] constexpr bool intersects(triangle<2, ScalarT> const& a, aabb<2, ScalarT> const& b)
@@ -854,11 +862,6 @@ template <class ScalarT>
 
     return true;
 }
-template <class ScalarT>
-[[nodiscard]] constexpr bool intersects(aabb<2, ScalarT> const& a, triangle<2, ScalarT> const& b)
-{
-    return intersects(b, a);
-}
 
 // NOTE: does NOT work for integer objects
 template <class ScalarT>
@@ -960,10 +963,5 @@ template <class ScalarT>
     // found no separating axis? -> intersection
     return true;
 }
-template <class ScalarT>
-[[nodiscard]] constexpr bool intersects(aabb<3, ScalarT> const& a, triangle<3, ScalarT> const& b)
-{
-    return intersects(b, a);
-}
 
 } // namespace tg
-- 
GitLab


From 93b920c4117cc40a0d5f9518748b0ba3284f5a91 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Sun, 8 Nov 2020 23:33:15 +0100
Subject: [PATCH 10/42] Added ray - halfspace intersection

---
 .../functions/objects/intersection.hh          | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index c8e0544..c0cb20d 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -214,6 +214,7 @@ constexpr optional<pos<D, ScalarT>> intersection(Obj const& obj, pos<D, ScalarT>
 template <class ScalarT>
 [[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<2, ScalarT> const& r, line<2, ScalarT> const& l)
 {
+    // r.origin + r.dir * t.x == l.pos + l.dir * t.y  <=>  (r.dir | -l.dir) * (t.x | t.y)^T == l.pos - r.origin
     auto M = mat<2, 2, ScalarT>::from_cols(r.dir, -l.dir);
     auto t = inverse(M) * (l.pos - r.origin);
     if (t.x < ScalarT(0))
@@ -267,6 +268,23 @@ template <int D, class ScalarT>
     return t;
 }
 
+// ray - halfspace
+template <int D, class ScalarT>
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, halfspace<D, ScalarT> const& h)
+{
+    // 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
+    auto dotND = dot(h.normal, r.dir);
+    if (dotND >= ScalarT(0))
+        return {};
+
+    return - dist / dotND;
+}
+
 // returns closest intersection point(s) of ray and sphere
 template <int D, class ScalarT>
 [[nodiscard]] constexpr optional<ScalarT> closest_intersection_parameter(ray<D, ScalarT> const& r, sphere<D, ScalarT> const& s)
-- 
GitLab


From 4b859ba49ce3b1f3667aa541563bec55e37c54cf Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Mon, 9 Nov 2020 13:48:22 +0100
Subject: [PATCH 11/42] Added ray2 - triangle2 intersection

---
 .../functions/objects/intersection.hh         | 29 ++++++++++++++++++-
 1 file changed, 28 insertions(+), 1 deletion(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index c0cb20d..4411f81 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -452,7 +452,34 @@ template <class ScalarT>
     return {hits, hit_cnt};
 }
 
-// ray - triangle
+// ray - triangle2
+template <class ScalarT>
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<2, ScalarT> const& r, triangle<2, ScalarT> const& t)
+{
+    ScalarT closestIntersection = tg::max<ScalarT>();
+    auto numIntersections = 0;
+    for (const auto& edge : edges_of(t))
+    {
+        const auto inter = intersection_parameter(r, edge);
+        if (inter.any())
+        {
+            numIntersections++;
+            closestIntersection = min(closestIntersection, inter.first());
+        }
+    }
+
+    if (numIntersections == 0)
+        return {}; // No intersection
+    if (numIntersections == 1)
+    {
+        // ray started within the triangle
+        TG_ASSERT(contains(t, r.origin));
+        return ScalarT(0);
+    }
+    return closestIntersection;
+}
+
+// ray - triangle3
 template <class ScalarT>
 [[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<3, ScalarT> const& r,
                                                                  triangle<3, ScalarT> const& t,
-- 
GitLab


From e7bfa13fcb36690637a923a9d4928898eab97b7c Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Mon, 9 Nov 2020 23:44:58 +0100
Subject: [PATCH 12/42] Fixed contains(halfspace3) shortcut selection

---
 src/typed-geometry/functions/objects/contains.hh | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/typed-geometry/functions/objects/contains.hh b/src/typed-geometry/functions/objects/contains.hh
index 123d9ef..46c1878 100644
--- a/src/typed-geometry/functions/objects/contains.hh
+++ b/src/typed-geometry/functions/objects/contains.hh
@@ -260,7 +260,7 @@ template <int D, class ScalarT>
 }
 
 template <int D, class ScalarT>
-[[nodiscard]] constexpr bool contains(halfspace<D, ScalarT> const& h, pos<2, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0))
+[[nodiscard]] constexpr bool contains(halfspace<D, ScalarT> const& h, pos<D, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0))
 {
     return signed_distance(p, h) <= eps;
 }
-- 
GitLab


From f80f3eeb06fa32865e283c080640f5d450062f7e Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Mon, 9 Nov 2020 23:55:14 +0100
Subject: [PATCH 13/42] Introduced ray_interval return type for solid
 intersections; Added ray aabb intersection

---
 .../functions/objects/intersection.hh         | 134 ++++++++++++++----
 src/typed-geometry/types/objects/aabb.hh      |   1 +
 2 files changed, 107 insertions(+), 28 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 4411f81..6bb4693 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -27,7 +27,7 @@
 // intersects(a, b)              -> bool
 // intersection(a, b)            -> ???
 // intersection_safe(a, b)       -> optional<???>
-// intersection_parameter(a, b)  -> coords? (for a ray: ray_hits<N, ScalarT>)
+// intersection_parameter(a, b)  -> coords? (for a ray: ray_hits<N, ScalarT> or optional<ray_interval> (when b is solid))
 // intersection_parameters(a, b) -> pair<coords, coords>?
 // intersection_exact(a, b)      -> variant
 // closest_intersection(a, b)            -> position/object (for a ray: optional<pos>)
@@ -48,6 +48,7 @@
 
 // Implementation guidelines:
 // explicit intersection_parameter(ray, obj), which gives closest_intersection_parameter and intersection position(s)
+//          intersection_parameter(ray, obj_boundary) already gives intersection_parameter(ray, obj) when boundary_of(obj) is defined and object is convex
 // explicit closest_intersection_parameter(ray, obj) if this is faster than computing all intersections
 // explicit intersects(obj, aabb), which gives intersects(aabb, obj)
 
@@ -103,6 +104,16 @@ private:
     HitT _hit[MaxHits];
 };
 
+/// describes a continuous interval on a ray between start and end
+template <class ScalarT>
+struct ray_interval
+{
+    ScalarT start;
+    ScalarT end;
+
+    [[nodiscard]] constexpr bool is_unbounded() const { return end == tg::max<ScalarT>(); }
+};
+
 
 // ====================================== Default Implementations ======================================
 // TODO: intersection_parameter from intersection_parameters
@@ -114,11 +125,11 @@ template <class A, class B>
     return intersection(a, b).has_value();
 }
 
-// if ray_hits intersection parameter is available and applicable, use that
+// if closest intersection parameter is available and applicable, use that
 template <int D, class ScalarT, class Obj>
-[[nodiscard]] constexpr auto intersects(ray<D, ScalarT> const& r, Obj const& obj) -> decltype(intersection_parameter(r, obj).any())
+[[nodiscard]] constexpr auto intersects(ray<D, ScalarT> const& r, Obj const& obj) -> decltype(closest_intersection_parameter(r, obj).has_value())
 {
-    return intersection_parameter(r, obj).any();
+    return closest_intersection_parameter(r, obj).has_value();
 }
 
 // parameters for intersects with aabb can switch order
@@ -178,6 +189,17 @@ template <int D, class ScalarT, class Obj>
     return {hits, ts.size()};
 }
 
+// if an optional ray_interval intersection parameter is available, use that
+template <int D, class ScalarT, class Obj, typename = std::enable_if_t<std::is_same_v<decltype(
+    intersection_parameter(std::declval<ray<D, ScalarT>>(), std::declval<Obj>())), optional<ray_interval<ScalarT>>>>>
+[[nodiscard]] constexpr optional<segment<D, ScalarT>> intersection(ray<D, ScalarT> const& r, Obj const& obj)
+{
+    auto ts = intersection_parameter(r, obj);
+    if (ts.has_value())
+        return {r[ts.value().start], r[ts.value().end]};
+    return {};
+}
+
 // if ray_hits intersection parameter is available, use that
 template <int D, class ScalarT, class Obj>
 [[nodiscard]] constexpr auto closest_intersection_parameter(ray<D, ScalarT> const& r, Obj const& obj)
@@ -189,6 +211,37 @@ template <int D, class ScalarT, class Obj>
     return {};
 }
 
+// if ray_interval intersection parameter is available, use that
+template <int D, class ScalarT, class Obj>
+[[nodiscard]] constexpr auto closest_intersection_parameter(ray<D, ScalarT> const& r, Obj const& obj)
+    -> enable_if<std::is_same_v<decltype(intersection_parameter(r, obj)), optional<ray_interval<ScalarT>>>, optional<ScalarT>>
+{
+    auto hits = intersection_parameter(r, obj);
+    if (hits.has_value())
+        return hits.value().start;
+    return {};
+}
+
+// if boundary_of a solid object returns ray_hits, use this to construct the ray_interval result of the solid intersection
+template <int D, class ScalarT, class Obj>
+[[nodiscard]] constexpr auto intersection_parameter(ray<D, ScalarT> const& r, Obj const& obj)
+    -> enable_if<std::is_same_v<decltype(intersection_parameter(r, boundary_of(obj))), ray_hits<2, ScalarT>>, optional<ray_interval<ScalarT>>>
+{
+    const auto inter = intersection_parameter(r, boundary_of(obj));
+
+    if (inter.size() == 0)
+        return {};
+
+    if (inter.size() == 1)
+    {
+        TG_ASSERT(contains(obj, r.origin));
+        return {{typename Obj::scalar_t(0), inter.first()}};
+    }
+
+    TG_ASSERT(inter.size() == 2);
+    return {{inter[0], inter[1]}};
+}
+
 // 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)
@@ -248,8 +301,8 @@ template <class ScalarT>
 template <int D, class ScalarT>
 [[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, plane<D, ScalarT> const& p)
 {
-    // if plane normal and ray direction are parallel there is no intersection
-    auto dotND = dot(p.normal, r.dir);
+    // if plane normal and ray direction are orthogonal there is no intersection
+    const auto dotND = dot(p.normal, r.dir);
     if (dotND == ScalarT(0))
         return {};
 
@@ -259,7 +312,7 @@ template <int D, class ScalarT>
     // <r.origin + t * r.dir, p.normal> = p.dis
     // t = (p.dis - <r.origin, p.normal>) / <r.dir, p.normal>
 
-    auto t = (p.dis - dot(p.normal, r.origin)) / dotND;
+    const auto t = (p.dis - dot(p.normal, r.origin)) / dotND;
 
     // check whether plane lies behind ray
     if (t < ScalarT(0))
@@ -270,7 +323,26 @@ template <int D, class ScalarT>
 
 // ray - halfspace
 template <int D, class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, halfspace<D, ScalarT> const& h)
+[[nodiscard]] constexpr optional<ray_interval<ScalarT>> intersection_parameter(ray<D, ScalarT> const& r, halfspace<D, ScalarT> const& h)
+{
+    // check if ray origin is already contained in the halfspace
+    const auto dotND = dot(h.normal, r.dir);
+    const auto dist = signed_distance(r.origin, h);
+    if (dist <= ScalarT(0))
+    {
+        if (dotND <= ScalarT(0)) // if ray point away from the surface, the entire ray intersects
+            return {{ScalarT(0), tg::max<ScalarT>()}};
+        return {{ScalarT(0), -dist / dotND}}; // otherwise the segment from ray origin until surface intersection is contained
+    }
+
+    // if ray points away from the halfspace there is no intersection
+    if (dotND >= ScalarT(0))
+        return {};
+
+    return {{-dist / dotND, tg::max<ScalarT>()}};
+}
+template <int D, class ScalarT>
+[[nodiscard]] constexpr optional<ScalarT> closest_intersection_parameter(ray<D, ScalarT> const& r, halfspace<D, ScalarT> const& h)
 {
     // check if ray origin is already contained in the halfspace
     const auto dist = signed_distance(r.origin, h);
@@ -278,37 +350,41 @@ template <int D, class ScalarT>
         return ScalarT(0);
 
     // if ray points away from the halfspace there is no intersection
-    auto dotND = dot(h.normal, r.dir);
+    const auto dotND = dot(h.normal, r.dir);
     if (dotND >= ScalarT(0))
         return {};
 
-    return - dist / dotND;
+    return -dist / dotND;
 }
 
-// returns closest intersection point(s) of ray and sphere
+// ray - aabb
 template <int D, class ScalarT>
-[[nodiscard]] constexpr optional<ScalarT> closest_intersection_parameter(ray<D, ScalarT> const& r, sphere<D, ScalarT> const& s)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, aabb_boundary<D, ScalarT> const& b)
 {
-    auto t = dot(s.center - r.origin, r.dir);
+    const auto tMin = (b.min - r.origin) / comp(r.dir);
+    const auto tMax = (b.max - r.origin) / comp(r.dir);
+    auto tFirst = max(min(tMin.x, tMax.x), min(tMin.y, tMax.y));
+    auto tSecond = min(max(tMin.x, tMax.x), max(tMin.y, tMax.y));
 
-    auto d_sqr = distance_sqr(r[t], s.center);
-    auto r_sqr = s.radius * s.radius;
-    if (d_sqr > r_sqr)
-        return {};
+    if constexpr (D >= 3)
+    {
+        tFirst = max(tFirst, min(tMin.z, tMax.z));
+        tSecond = min(tSecond, max(tMin.z, tMax.z));
+    }
 
-    auto dt = sqrt(r_sqr - d_sqr);
+    // No intersection if ray starts behind aabb or if it misses the aabb
+    if (tSecond < ScalarT(0) || tFirst > tSecond)
+        return {};
 
-    if (t - dt >= 0)
-        return t - dt;
-    if (t + dt >= 0)
-        return t + dt;
+    if (tFirst < ScalarT(0))
+        return tSecond;
 
-    return {};
+    return {tFirst, tSecond};
 }
 
 // returns intersection point(s) of ray and sphere
 template <int D, class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, sphere<D, ScalarT> const& s)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, sphere_boundary<D, ScalarT> const& s)
 {
     auto t = dot(s.center - r.origin, r.dir);
 
@@ -454,9 +530,10 @@ template <class ScalarT>
 
 // ray - triangle2
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<2, ScalarT> const& r, triangle<2, ScalarT> const& t)
+[[nodiscard]] constexpr optional<ray_interval<ScalarT>> intersection_parameter(ray<2, ScalarT> const& r, triangle<2, ScalarT> const& t)
 {
     ScalarT closestIntersection = tg::max<ScalarT>();
+    ScalarT furtherIntersection = tg::min<ScalarT>();
     auto numIntersections = 0;
     for (const auto& edge : edges_of(t))
     {
@@ -465,6 +542,7 @@ template <class ScalarT>
         {
             numIntersections++;
             closestIntersection = min(closestIntersection, inter.first());
+            furtherIntersection = max(furtherIntersection, inter.first());
         }
     }
 
@@ -474,9 +552,9 @@ template <class ScalarT>
     {
         // ray started within the triangle
         TG_ASSERT(contains(t, r.origin));
-        return ScalarT(0);
+        return {{ScalarT(0), furtherIntersection}};
     }
-    return closestIntersection;
+    return {{closestIntersection, furtherIntersection}};
 }
 
 // ray - triangle3
@@ -517,7 +595,7 @@ template <class ScalarT>
 
 // ray - box
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, box<3, ScalarT> const& b)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, box_boundary<3, ScalarT> const& b)
 {
     // see https://github.com/gszauer/GamePhysicsCookbook/blob/master/Code/Geometry3D.cpp
 
diff --git a/src/typed-geometry/types/objects/aabb.hh b/src/typed-geometry/types/objects/aabb.hh
index b0e9121..5afc17d 100644
--- a/src/typed-geometry/types/objects/aabb.hh
+++ b/src/typed-geometry/types/objects/aabb.hh
@@ -48,6 +48,7 @@ using aabb_boundary = aabb<D, ScalarT, boundary_tag>;
 template <int D, class ScalarT, class TraitsT>
 struct aabb
 {
+    using scalar_t = ScalarT;
     using vec_t = vec<D, ScalarT>;
     using pos_t = pos<D, ScalarT>;
 
-- 
GitLab


From f0fd7928906fffa2e12742e2fc5254cd74b91107 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Tue, 10 Nov 2020 00:43:31 +0100
Subject: [PATCH 14/42] Improved ray - aabb intersection, now for all
 dimensions

---
 .../functions/objects/intersection.hh         | 23 +++++++++++--------
 1 file changed, 13 insertions(+), 10 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 6bb4693..ed8a6eb 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -4,8 +4,10 @@
 #include <typed-geometry/feature/assert.hh>
 #include <typed-geometry/functions/basic/scalar_math.hh>
 
-#include <typed-geometry/types/objects/capsule.hh>
+#include <typed-geometry/types/objects/aabb.hh>
 #include <typed-geometry/types/objects/cylinder.hh>
+#include <typed-geometry/types/objects/halfspace.hh>
+#include <typed-geometry/types/objects/line.hh>
 #include <typed-geometry/types/objects/plane.hh>
 #include <typed-geometry/types/objects/ray.hh>
 #include <typed-geometry/types/objects/segment.hh>
@@ -17,7 +19,6 @@
 #include <typed-geometry/functions/vector/length.hh>
 
 #include "aabb.hh"
-#include "closest_points.hh"
 #include "contains.hh"
 #include "direction.hh"
 #include "normal.hh"
@@ -361,15 +362,17 @@ template <int D, class ScalarT>
 template <int D, class ScalarT>
 [[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, aabb_boundary<D, ScalarT> const& b)
 {
-    const auto tMin = (b.min - r.origin) / comp(r.dir);
-    const auto tMax = (b.max - r.origin) / comp(r.dir);
-    auto tFirst = max(min(tMin.x, tMax.x), min(tMin.y, tMax.y));
-    auto tSecond = min(max(tMin.x, tMax.x), max(tMin.y, tMax.y));
-
-    if constexpr (D >= 3)
+    auto tFirst = tg::min<ScalarT>();
+    auto tSecond = tg::max<ScalarT>();
+    for (auto i = 0; i < D; ++i)
     {
-        tFirst = max(tFirst, min(tMin.z, tMax.z));
-        tSecond = min(tSecond, max(tMin.z, tMax.z));
+        if (abs(r.dir[i]) > ScalarT(100) * tg::epsilon<ScalarT>)
+        {
+            const auto tMin = (b.min[i] - r.origin[i]) / r.dir[i];
+            const auto tMax = (b.max[i] - r.origin[i]) / r.dir[i];
+            tFirst = max(tFirst, min(tMin, tMax));
+            tSecond = min(tSecond, max(tMin, tMax));
+        }
     }
 
     // No intersection if ray starts behind aabb or if it misses the aabb
-- 
GitLab


From f4de456f3306de8f769c989cf861a20fa7d742bf Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Tue, 10 Nov 2020 01:33:29 +0100
Subject: [PATCH 15/42] Fixed ray - aabb intersection for axis aligned rays

---
 src/typed-geometry/functions/objects/intersection.hh | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index ed8a6eb..14f8880 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -373,9 +373,11 @@ template <int D, class ScalarT>
             tFirst = max(tFirst, min(tMin, tMax));
             tSecond = min(tSecond, max(tMin, tMax));
         }
+        else if (r.origin[i] < b.min[i] || r.origin[i] > b.max[i])
+            return {}; // ray parallel to this axis and outside of aabb
     }
 
-    // No intersection if ray starts behind aabb or if it misses the aabb
+    // no intersection if ray starts behind aabb or if it misses the aabb
     if (tSecond < ScalarT(0) || tFirst > tSecond)
         return {};
 
-- 
GitLab


From 3a9f8048463c4ce8197fae53a12d2c4e0f151454 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Tue, 10 Nov 2020 03:09:30 +0100
Subject: [PATCH 16/42] Reimplemented ray - box intersection

---
 .../functions/objects/intersection.hh         | 95 +++++++------------
 1 file changed, 35 insertions(+), 60 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 14f8880..17d983c 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -5,6 +5,7 @@
 #include <typed-geometry/functions/basic/scalar_math.hh>
 
 #include <typed-geometry/types/objects/aabb.hh>
+#include <typed-geometry/types/objects/box.hh>
 #include <typed-geometry/types/objects/cylinder.hh>
 #include <typed-geometry/types/objects/halfspace.hh>
 #include <typed-geometry/types/objects/line.hh>
@@ -374,10 +375,42 @@ template <int D, class ScalarT>
             tSecond = min(tSecond, max(tMin, tMax));
         }
         else if (r.origin[i] < b.min[i] || r.origin[i] > b.max[i])
-            return {}; // ray parallel to this axis and outside of aabb
+            return {}; // ray parallel to this axis and outside of the aabb
     }
 
-    // no intersection if ray starts behind aabb or if it misses the aabb
+    // no intersection if ray starts behind the aabb or if it misses the aabb
+    if (tSecond < ScalarT(0) || tFirst > tSecond)
+        return {};
+
+    if (tFirst < ScalarT(0))
+        return tSecond;
+
+    return {tFirst, tSecond};
+}
+
+// ray - box
+template <int D, class ScalarT>
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, box_boundary<D, ScalarT> const& b)
+{
+    const auto bMin = b[comp<D, ScalarT>(-1)] - r.origin;
+    const auto bMax = b[comp<D, ScalarT>(1)] - r.origin;
+    auto tFirst = tg::min<ScalarT>();
+    auto tSecond = tg::max<ScalarT>();
+    for (auto i = 0; i < D; ++i)
+    {
+        const auto rDir = dot(r.dir, b.half_extents[i]);
+        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
+    }
+
+    // no intersection if ray starts behind the box or if it misses the box
     if (tSecond < ScalarT(0) || tFirst > tSecond)
         return {};
 
@@ -598,64 +631,6 @@ template <class ScalarT>
     return (lambda > ScalarT(0)) ? lambda : ray_hits<1, ScalarT>();
 }
 
-// ray - box
-template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, box_boundary<3, ScalarT> const& b)
-{
-    // see https://github.com/gszauer/GamePhysicsCookbook/blob/master/Code/Geometry3D.cpp
-
-    auto const p = b.center - r.origin;
-
-    // TODO: without normalize
-    auto const X = normalize(b.half_extents[0]);
-    auto const Y = normalize(b.half_extents[1]);
-    auto const Z = normalize(b.half_extents[2]);
-
-    auto f = tg::vec3(dot(X, r.dir), //
-                      dot(Y, r.dir), //
-                      dot(Z, r.dir));
-
-    auto const e = tg::vec3(dot(X, p), //
-                            dot(Y, p), //
-                            dot(Z, p));
-
-    auto const size = tg::vec3(length(b.half_extents[0]), //
-                               length(b.half_extents[1]), //
-                               length(b.half_extents[2]));
-
-    float t[6] = {0, 0, 0, 0, 0, 0};
-    for (int i = 0; i < 3; ++i)
-    {
-        if (abs(f[i]) < tg::epsilon<ScalarT> / 1000)
-        {
-            if (-e[i] - size[i] > 0 || -e[i] + size[i] < 0)
-                return {};
-
-            f[i] = tg::epsilon<ScalarT>; // Avoid div by 0!
-        }
-
-        t[i * 2 + 0] = (e[i] + size[i]) / f[i]; // tmin[x, y, z]
-        t[i * 2 + 1] = (e[i] - size[i]) / f[i]; // tmax[x, y, z]
-    }
-
-    auto const tmin = tg::max(tg::max(tg::min(t[0], t[1]), tg::min(t[2], t[3])), tg::min(t[4], t[5]));
-    auto const tmax = tg::min(tg::min(tg::max(t[0], t[1]), tg::max(t[2], t[3])), tg::max(t[4], t[5]));
-
-    // if tmax < 0, ray is intersecting AABB
-    // but entire AABB is behing it's origin
-    if (tmax < 0)
-        return {};
-
-    // if tmin > tmax, ray doesn't intersect AABB
-    if (tmin > tmax)
-        return {};
-
-    if (tmin >= 0) // two valid intersections
-        return {tmin, tmax};
-    else // one valid intersection
-        return tmax;
-}
-
 
 // ====================================== Object - Object Intersections ======================================
 
-- 
GitLab


From 12ebf69c1c78bd298a5abf9cd71204142bc5dfbd Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Wed, 11 Nov 2020 00:34:58 +0100
Subject: [PATCH 17/42] Added helper functions to ease implementation of
 intersections with compound objects

---
 .../functions/objects/intersection.hh         | 72 ++++++++++++++-----
 1 file changed, 53 insertions(+), 19 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 17d983c..ca1ae81 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -53,6 +53,10 @@
 //          intersection_parameter(ray, obj_boundary) already gives intersection_parameter(ray, obj) when boundary_of(obj) is defined and object is convex
 // explicit closest_intersection_parameter(ray, obj) if this is faster than computing all intersections
 // explicit intersects(obj, aabb), which gives intersects(aabb, obj)
+//
+// for convex compound objects (like cylinder or pyramids), decompose the object into primitive shapes and pass them to a helper function:
+// - call merge_hits(ray, objPart1, objPart2, ...) in the implementation of intersection_parameter(ray, obj_boundary)
+// - call intersects_any(ray, objPart1, objPart2, ...) in the implementation of intersects(ray, obj<TraitsT>), which can shortcut and be faster than the default
 
 
 namespace tg
@@ -263,6 +267,42 @@ constexpr optional<pos<D, ScalarT>> intersection(Obj const& obj, pos<D, ScalarT>
 }
 
 
+// ====================================== Helper functions ======================================
+
+namespace detail
+{
+// intersects the given ray with all given objects and returns the concatenated intersections. A maximal number of 2 intersections is assumed.
+template <int D, class ScalarT, class... ObjsT>
+[[nodiscard]] constexpr ray_hits<2, ScalarT> merge_hits(ray<D, ScalarT> const& ray, ObjsT const&... objs)
+{
+    ScalarT hits[2];
+    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(ray, 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};
+}
+
+// returns true, iff the given ray intersects any of the given objects (with short-circuiting after the first intersection)
+template <int D, class ScalarT, class... ObjsT>
+[[nodiscard]] constexpr bool intersects_any(ray<D, ScalarT> const& ray, ObjsT const&... objs)
+{
+    return (intersects(ray, objs) || ...);
+}
+}
+
 // ====================================== Ray - Object Intersections ======================================
 
 // ray - line
@@ -493,27 +533,21 @@ template <class ScalarT>
 
 // ray - cylinder
 template <class ScalarT>
-[[nodiscard]] constexpr optional<ScalarT> closest_intersection_parameter(ray<3, ScalarT> const& r, cylinder_boundary<3, ScalarT> const& c)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, cylinder_boundary<3, ScalarT> const& c)
 {
     const auto caps = caps_of(c);
-    const auto t_cap0 = intersection_parameter(r, caps[0]);
-    const auto t_cap1 = intersection_parameter(r, caps[1]);
-    const auto t_cyl = closest_intersection_parameter(r, boundary_no_caps_of(c));
-
-    optional<ScalarT> t;
-
-    if (t_cyl.has_value())
-        t = t_cyl.value();
-
-    if (t_cap0.any() && (!t.has_value() || t.value() > t_cap0.first()))
-        t = t_cap0.first();
-
-    if (t_cap1.any() && (!t.has_value() || t.value() > t_cap1.first()))
-        t = t_cap1.first();
-
-    TG_INTERNAL_ASSERT(!t.has_value() || t.value() >= 0);
-
-    return t;
+    return detail::merge_hits(r, caps[0], caps[1], boundary_no_caps_of(c));
+}
+template <class ScalarT, class TraitsT>
+[[nodiscard]] constexpr bool intersects(ray<3, ScalarT> const& r, cylinder<3, ScalarT, TraitsT> const& c)
+{
+    if constexpr (std::is_same_v<TraitsT, boundary_no_caps_tag>)
+        return intersection_parameter(r, c).any();
+    else
+    {
+        const auto caps = caps_of(c);
+        return detail::intersects_any(r, boundary_no_caps_of(c), caps[0], caps[1]);
+    }
 }
 
 // ray - tube
-- 
GitLab


From 938952c7b81af3e7f820a24e05edd5e46448a0aa Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Wed, 11 Nov 2020 00:35:52 +0100
Subject: [PATCH 18/42] Added contains(capsule), which is faster than the
 default implementation

---
 src/typed-geometry/functions/objects/contains.hh | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/src/typed-geometry/functions/objects/contains.hh b/src/typed-geometry/functions/objects/contains.hh
index 46c1878..25cecc8 100644
--- a/src/typed-geometry/functions/objects/contains.hh
+++ b/src/typed-geometry/functions/objects/contains.hh
@@ -4,6 +4,7 @@
 
 #include <typed-geometry/types/objects/aabb.hh>
 #include <typed-geometry/types/objects/box.hh>
+#include <typed-geometry/types/objects/capsule.hh>
 #include <typed-geometry/types/objects/cone.hh>
 #include <typed-geometry/types/objects/cylinder.hh>
 #include <typed-geometry/types/objects/ellipse.hh>
@@ -221,6 +222,19 @@ template <class ScalarT, class TraitsT>
     return x2 / pow2(a + eps) + y2 / pow2(b + eps) + z2 / pow2(c + eps) <= ScalarT(1);
 }
 
+template <class ScalarT>
+[[nodiscard]] constexpr bool contains(capsule<3, ScalarT> const& c, pos<3, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0))
+{
+    auto r = c.radius + eps;
+    return distance_sqr(c.axis, p) <= pow2(r);
+}
+template <class ScalarT>
+[[nodiscard]] constexpr bool contains(capsule_boundary<3, ScalarT> const& c, pos<3, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0))
+{
+    auto d2 = distance_sqr(c.axis, p);
+    return pow2(max(ScalarT(0), c.radius - eps)) <= d2 && d2 <= pow2(c.radius + eps);
+}
+
 template <int D, class ScalarT>
 [[nodiscard]] constexpr bool contains(sphere<D, ScalarT> const& s, pos<D, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0))
 {
-- 
GitLab


From 834104049d6de589530b0c4fc05c33279b6a9a05 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Wed, 11 Nov 2020 02:11:40 +0100
Subject: [PATCH 19/42] Implemented ray intersections for hemisphere and
 capsule

---
 .../functions/objects/intersection.hh         | 60 ++++++++++++++++++-
 1 file changed, 59 insertions(+), 1 deletion(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index ca1ae81..20e8c55 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -6,8 +6,10 @@
 
 #include <typed-geometry/types/objects/aabb.hh>
 #include <typed-geometry/types/objects/box.hh>
+#include <typed-geometry/types/objects/capsule.hh>
 #include <typed-geometry/types/objects/cylinder.hh>
 #include <typed-geometry/types/objects/halfspace.hh>
+#include <typed-geometry/types/objects/hemisphere.hh>
 #include <typed-geometry/types/objects/line.hh>
 #include <typed-geometry/types/objects/plane.hh>
 #include <typed-geometry/types/objects/ray.hh>
@@ -21,6 +23,7 @@
 
 #include "aabb.hh"
 #include "contains.hh"
+#include "coordinates.hh"
 #include "direction.hh"
 #include "normal.hh"
 
@@ -305,6 +308,16 @@ template <int D, class ScalarT, class... ObjsT>
 
 // ====================================== Ray - Object Intersections ======================================
 
+// ray - point
+template <class ScalarT>
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<1, ScalarT> const& r, pos<1, ScalarT> const& p)
+{
+    const auto t = coordinates(r, p);
+    if (t >= ScalarT(0))
+        return t;
+    return {};
+}
+
 // ray - line
 template <class ScalarT>
 [[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<2, ScalarT> const& r, line<2, ScalarT> const& l)
@@ -497,6 +510,34 @@ template <class ScalarT>
     return t;
 }
 
+// ray - hemisphere
+template <int D, class ScalarT>
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, hemisphere_boundary<D, ScalarT> const& h)
+{
+    return detail::merge_hits(r, caps_of(h), boundary_no_caps_of(h));
+}
+template <int D, class ScalarT>
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, hemisphere_boundary_no_caps<D, ScalarT> const& h)
+{
+    ScalarT hits[2];
+    auto numHits = 0;
+    const auto sphereHits = intersection_parameter(r, sphere_boundary<D, ScalarT>(h.center, h.radius));
+    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)
+        if (contains(halfSpace, r[hit]))
+            hits[numHits++] = hit;
+
+    return {hits, numHits};
+}
+template <int D, class ScalarT, class TraitsT>
+[[nodiscard]] constexpr bool intersects(ray<D, ScalarT> const& r, hemisphere<D, ScalarT, TraitsT> const& h)
+{
+    if constexpr (std::is_same_v<TraitsT, boundary_no_caps_tag>)
+        return intersection_parameter(r, h).any();
+    else
+        return detail::intersects_any(r, caps_of(h), boundary_no_caps_of(h));
+}
+
 // ray - quadric (as an isosurface, not error quadric)
 template <class ScalarT>
 [[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, quadric<3, ScalarT> const& q)
@@ -531,6 +572,23 @@ template <class ScalarT>
     return {};
 }
 
+// ray - capsule
+template <class ScalarT>
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, capsule_boundary<3, ScalarT> const& c)
+{
+    using caps_t = hemisphere_boundary_no_caps<3, ScalarT>;
+    const auto n = direction(c);
+    return detail::merge_hits(r, caps_t(c.axis.pos0, c.radius, -n), caps_t(c.axis.pos1, c.radius, n),
+                              cylinder_boundary_no_caps<3, ScalarT>(c.axis, c.radius));
+}
+template <class ScalarT, class TraitsT>
+[[nodiscard]] constexpr bool intersects(ray<3, ScalarT> const& r, capsule<3, ScalarT, TraitsT> const& c)
+{
+    using caps_t = sphere_boundary<3, ScalarT>; // spheres are faster than hemispheres and equivalent for the yes/no decision
+    return detail::intersects_any(r, caps_t(c.axis.pos0, c.radius), caps_t(c.axis.pos1, c.radius),
+                                  cylinder_boundary_no_caps<3, ScalarT>(c.axis, c.radius));
+}
+
 // ray - cylinder
 template <class ScalarT>
 [[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, cylinder_boundary<3, ScalarT> const& c)
@@ -546,7 +604,7 @@ template <class ScalarT, class TraitsT>
     else
     {
         const auto caps = caps_of(c);
-        return detail::intersects_any(r, boundary_no_caps_of(c), caps[0], caps[1]);
+        return detail::intersects_any(r, caps[0], caps[1], boundary_no_caps_of(c));
     }
 }
 
-- 
GitLab


From 5153d9f86bac6054ab2097aa7183ada7ee891a03 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Wed, 11 Nov 2020 02:19:44 +0100
Subject: [PATCH 20/42] Moved ray - quadric intersection code

---
 .../functions/objects/intersection.hh         | 68 +++++++++----------
 1 file changed, 34 insertions(+), 34 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 20e8c55..04a6d61 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -538,40 +538,6 @@ template <int D, class ScalarT, class TraitsT>
         return detail::intersects_any(r, caps_of(h), boundary_no_caps_of(h));
 }
 
-// ray - quadric (as an isosurface, not error quadric)
-template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, quadric<3, ScalarT> const& q)
-{
-    const auto Ad = q.A() * r.dir;
-    const auto p = r.origin;
-
-    // Substituting x in Quadric equation x^TAx + 2b^Tx + c = 0 by ray equation x = t * dir + p yields
-    // d^TAd t^2 + (2p^TAd + 2bd) t + p^TAp + 2bp + c = 0
-    const auto a = dot(r.dir, Ad);
-    const auto b = ScalarT(2) * (dot(p, Ad) + dot(q.b(), r.dir));
-    const auto c = dot(p, q.A() * vec3(p)) + ScalarT(2) * dot(q.b(), p) + q.c;
-
-    // Solve the quadratic equation ax^2 + bx + c = 0
-    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);
-
-    auto tMin = min(t1, t2);
-    auto tMax = max(t1, t2);
-
-    if (tMin >= ScalarT(0))
-        return {tMin, tMax};
-
-    if (tMax >= ScalarT(0))
-        return tMax;
-
-    return {};
-}
-
 // ray - capsule
 template <class ScalarT>
 [[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, capsule_boundary<3, ScalarT> const& c)
@@ -723,6 +689,40 @@ template <class ScalarT>
     return (lambda > ScalarT(0)) ? lambda : ray_hits<1, ScalarT>();
 }
 
+// ray - quadric_boundary (as an isosurface, not error quadric)
+template <class ScalarT>
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, quadric<3, ScalarT> const& q)
+{
+    const auto Ad = q.A() * r.dir;
+    const auto p = r.origin;
+
+    // Substituting x in Quadric equation x^TAx + 2b^Tx + c = 0 by ray equation x = t * dir + p yields
+    // d^TAd t^2 + (2p^TAd + 2bd) t + p^TAp + 2bp + c = 0
+    const auto a = dot(r.dir, Ad);
+    const auto b = ScalarT(2) * (dot(p, Ad) + dot(q.b(), r.dir));
+    const auto c = dot(p, q.A() * vec3(p)) + ScalarT(2) * dot(q.b(), p) + q.c;
+
+    // Solve the quadratic equation ax^2 + bx + c = 0
+    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);
+
+    auto tMin = min(t1, t2);
+    auto tMax = max(t1, t2);
+
+    if (tMin >= ScalarT(0))
+        return {tMin, tMax};
+
+    if (tMax >= ScalarT(0))
+        return tMax;
+
+    return {};
+}
+
 
 // ====================================== Object - Object Intersections ======================================
 
-- 
GitLab


From 554dc0264512b88aaa64c3215bc3bcd9d8f90492 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Mon, 16 Nov 2020 10:45:29 +0100
Subject: [PATCH 21/42] Trying to fix template instantiation

---
 src/typed-geometry/functions/objects/intersection.hh | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 04a6d61..64cb693 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -53,7 +53,8 @@
 
 // Implementation guidelines:
 // explicit intersection_parameter(ray, obj), which gives closest_intersection_parameter and intersection position(s)
-//          intersection_parameter(ray, obj_boundary) already gives intersection_parameter(ray, obj) when boundary_of(obj) is defined and object is convex
+//          intersection_parameter(ray, obj_boundary) already gives intersection_parameter(ray, obj) when boundary_of(obj) is defined and object is convex and finite
+//          for infinite objects, a separate solid version is needed when the boundary version does not distinguish between completely contained and completely outside
 // explicit closest_intersection_parameter(ray, obj) if this is faster than computing all intersections
 // explicit intersects(obj, aabb), which gives intersects(aabb, obj)
 //
@@ -234,7 +235,7 @@ template <int D, class ScalarT, class Obj>
 // if boundary_of a solid object returns ray_hits, use this to construct the ray_interval result of the solid intersection
 template <int D, class ScalarT, class Obj>
 [[nodiscard]] constexpr auto intersection_parameter(ray<D, ScalarT> const& r, Obj const& obj)
-    -> enable_if<std::is_same_v<decltype(intersection_parameter(r, boundary_of(obj))), ray_hits<2, ScalarT>>, optional<ray_interval<ScalarT>>>
+    -> enable_if<std::conjunction_v<std::negation<std::is_same<Obj, decltype(boundary_of(obj))>>, std::is_same<decltype(intersection_parameter(r, boundary_of(obj))), ray_hits<2, ScalarT>>>, optional<ray_interval<ScalarT>>>
 {
     const auto inter = intersection_parameter(r, boundary_of(obj));
 
-- 
GitLab


From 39b80b3765f73c79a5673a6a1e01fd62f4bc34a4 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Mon, 16 Nov 2020 12:32:27 +0100
Subject: [PATCH 22/42] Fixed template recursion

---
 src/typed-geometry/functions/objects/intersection.hh | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 64cb693..7440307 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -235,7 +235,7 @@ template <int D, class ScalarT, class Obj>
 // if boundary_of a solid object returns ray_hits, use this to construct the ray_interval result of the solid intersection
 template <int D, class ScalarT, class Obj>
 [[nodiscard]] constexpr auto intersection_parameter(ray<D, ScalarT> const& r, Obj const& obj)
-    -> enable_if<std::conjunction_v<std::negation<std::is_same<Obj, decltype(boundary_of(obj))>>, std::is_same<decltype(intersection_parameter(r, boundary_of(obj))), ray_hits<2, ScalarT>>>, optional<ray_interval<ScalarT>>>
+    -> enable_if<!std::is_same_v<Obj, decltype(boundary_of(obj))>, optional<ray_interval<ScalarT>>>
 {
     const auto inter = intersection_parameter(r, boundary_of(obj));
 
-- 
GitLab


From d0f7c669e29e5d8a2de973c614563eb9a352a6ee Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Mon, 16 Nov 2020 12:56:24 +0100
Subject: [PATCH 23/42] Fixed initializer list

---
 src/typed-geometry/functions/objects/intersection.hh | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 7440307..d6b08aa 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -206,7 +206,7 @@ template <int D, class ScalarT, class Obj, typename = std::enable_if_t<std::is_s
 {
     auto ts = intersection_parameter(r, obj);
     if (ts.has_value())
-        return {r[ts.value().start], r[ts.value().end]};
+        return {{r[ts.value().start], r[ts.value().end]}};
     return {};
 }
 
@@ -237,7 +237,7 @@ template <int D, class ScalarT, class Obj>
 [[nodiscard]] constexpr auto intersection_parameter(ray<D, ScalarT> const& r, Obj const& obj)
     -> enable_if<!std::is_same_v<Obj, decltype(boundary_of(obj))>, optional<ray_interval<ScalarT>>>
 {
-    const auto inter = intersection_parameter(r, boundary_of(obj));
+    const ray_hits<2, ScalarT> inter = intersection_parameter(r, boundary_of(obj));
 
     if (inter.size() == 0)
         return {};
-- 
GitLab


From 1a09d9b1c877e111a044535ff42228dcdd6f1f73 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Mon, 16 Nov 2020 13:30:28 +0100
Subject: [PATCH 24/42] Test for another CI problem

---
 src/typed-geometry/detail/optional.hh | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/typed-geometry/detail/optional.hh b/src/typed-geometry/detail/optional.hh
index f7f3fd7..b8b8965 100644
--- a/src/typed-geometry/detail/optional.hh
+++ b/src/typed-geometry/detail/optional.hh
@@ -47,6 +47,9 @@ public:
 template <class T, class Stream>
 auto operator<<(Stream& os, optional<T> const& o) -> decltype(os << o.value(), os << "", os)
 {
-    return o.has_value() ? (os << o.value()) : (os << "[empty]");
+    if (o.has_value())
+        return (os << o.value());
+    else
+        return (os << "[empty]");
 }
 }
-- 
GitLab


From 148405e7bfc301ae3fe5b7fe7e65edbde956c3f8 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Mon, 16 Nov 2020 14:10:36 +0100
Subject: [PATCH 25/42] Workaround for gcc bug

---
 src/typed-geometry/detail/optional.hh | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/src/typed-geometry/detail/optional.hh b/src/typed-geometry/detail/optional.hh
index b8b8965..dba24b1 100644
--- a/src/typed-geometry/detail/optional.hh
+++ b/src/typed-geometry/detail/optional.hh
@@ -45,11 +45,8 @@ public:
 };
 
 template <class T, class Stream>
-auto operator<<(Stream& os, optional<T> const& o) -> decltype(os << o.value(), os << "", os)
+auto operator<<(Stream& os, optional<T> const& o) -> decltype(os << o.value(), os << std::declval<char const*>(), os)
 {
-    if (o.has_value())
-        return (os << o.value());
-    else
-        return (os << "[empty]");
+    return o.has_value() ? (os << o.value()) : (os << "[empty]");
 }
 }
-- 
GitLab


From 207fc1b1aaae8a6e63f4c0557b17e0ea974b262f Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Mon, 16 Nov 2020 15:27:29 +0100
Subject: [PATCH 26/42] Added ray - box2in3 intersection

---
 .../functions/objects/intersection.hh         | 45 ++++++++++++-------
 1 file changed, 29 insertions(+), 16 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index d6b08aa..fdd0163 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -473,8 +473,36 @@ template <int D, class ScalarT>
 
     return {tFirst, tSecond};
 }
+template <class ScalarT>
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, box<2, ScalarT, 3> const& b)
+{
+    const auto t = intersection_parameter(r, plane<3, ScalarT>(normal_of(b), b.center));
+    if (!t.any())
+        return {};
+
+    const auto p = r[t.first()] - b.center;
+    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;
+}
 
-// returns intersection point(s) of ray and sphere
+// ray - disk
+template <class ScalarT>
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, sphere<2, ScalarT, 3> const& d)
+{
+    const auto t = intersection_parameter(r, plane<3, ScalarT>(d.normal, d.center));
+    if (!t.any())
+        return {};
+
+    const auto p = r[t.first()];
+    if (distance_sqr(p, d.center) > d.radius * d.radius)
+        return {};
+
+    return t;
+}
+
+// ray - sphere
 template <int D, class ScalarT>
 [[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, sphere_boundary<D, ScalarT> const& s)
 {
@@ -496,21 +524,6 @@ template <int D, class ScalarT>
     return {};
 }
 
-// ray - disk
-template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, sphere<2, ScalarT, 3> const& d)
-{
-    auto const t = intersection_parameter(r, plane<3, ScalarT>(d.normal, d.center));
-    if (!t.any())
-        return {};
-
-    auto const p = r[t.first()];
-    if (distance_sqr(p, d.center) > d.radius * d.radius)
-        return {};
-
-    return t;
-}
-
 // ray - hemisphere
 template <int D, class ScalarT>
 [[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, hemisphere_boundary<D, ScalarT> const& h)
-- 
GitLab


From 85411dd35f064a74ec65a84bd630e9a2c07e6f0e Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Mon, 16 Nov 2020 15:28:54 +0100
Subject: [PATCH 27/42] Fixed inf_cylinder description; Added inf_tube
 shorthands

---
 src/typed-geometry/types/objects/inf_cylinder.hh |  8 ++++----
 src/typed-geometry/types/objects/tube.hh         | 12 ++++++++++++
 2 files changed, 16 insertions(+), 4 deletions(-)

diff --git a/src/typed-geometry/types/objects/inf_cylinder.hh b/src/typed-geometry/types/objects/inf_cylinder.hh
index 8d1aa96..2a97962 100644
--- a/src/typed-geometry/types/objects/inf_cylinder.hh
+++ b/src/typed-geometry/types/objects/inf_cylinder.hh
@@ -8,16 +8,16 @@
 namespace tg
 {
 /**
- * infinite tube
+ * infinite cylinder
  *
- * An infinite tube is defined by a position and an axis
+ * An infinite cylinder is defined by an axis and a radius
  */
 template <int D, class ScalarT, class TraitsT = default_object_tag>
 struct inf_cylinder;
 
-// Common infinte tube types
+// Common infinite cylinder types
 using inf_cylinder3 = inf_cylinder<3, f32>;
-using finttube3 = inf_cylinder<3, f32>;
+using finf_cylinder3 = inf_cylinder<3, f32>;
 using dinf_cylinder3 = inf_cylinder<3, f64>;
 using iinf_cylinder3 = inf_cylinder<3, i32>;
 using uinf_cylinder3 = inf_cylinder<3, u32>;
diff --git a/src/typed-geometry/types/objects/tube.hh b/src/typed-geometry/types/objects/tube.hh
index bc90acd..16e1f54 100644
--- a/src/typed-geometry/types/objects/tube.hh
+++ b/src/typed-geometry/types/objects/tube.hh
@@ -1,6 +1,7 @@
 #pragma once
 
 #include "cylinder.hh"
+#include "inf_cylinder.hh"
 
 namespace tg
 {
@@ -13,4 +14,15 @@ using ftube3 = tube<3, f32>;
 using dtube3 = tube<3, f64>;
 using itube3 = tube<3, i32>;
 using utube3 = tube<3, u32>;
+
+
+template <int D, class ScalarT>
+using inf_tube = inf_cylinder_boundary<D, ScalarT>;
+
+// Common inf_tube types
+using inf_tube3 = inf_tube<3, f32>;
+using finf_tube3 = inf_tube<3, f32>;
+using dinf_tube3 = inf_tube<3, f64>;
+using iinf_tube3 = inf_tube<3, i32>;
+using uinf_tube3 = inf_tube<3, u32>;
 }
-- 
GitLab


From 8a4c563907bb83df99054881f48002c789472629 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Mon, 16 Nov 2020 23:48:00 +0100
Subject: [PATCH 28/42] Added minmax; Added ray - inf_cylinder2_boundary
 intersection

---
 .../functions/basic/scalar_math.hh            |  8 ++++++
 .../functions/objects/intersection.hh         | 26 ++++++++++++++++---
 .../types/objects/inf_cylinder.hh             |  1 +
 3 files changed, 31 insertions(+), 4 deletions(-)

diff --git a/src/typed-geometry/functions/basic/scalar_math.hh b/src/typed-geometry/functions/basic/scalar_math.hh
index 59f4d24..ef841c5 100644
--- a/src/typed-geometry/functions/basic/scalar_math.hh
+++ b/src/typed-geometry/functions/basic/scalar_math.hh
@@ -144,6 +144,14 @@ template <class A, class B>
     return a < b ? b : a;
 }
 
+template <class T>
+[[nodiscard]] constexpr pair<T, T> minmax(T const& a, T const& b)
+{
+    if (b < a)
+        return {b, a};
+    return {a, b};
+}
+
 template <class A, class B, class C>
 [[nodiscard]] constexpr auto clamp(A const& a, B const& min_value, C const& max_value) -> decltype(min(max(a, min_value), max_value))
 {
diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index fdd0163..dd25a67 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -587,8 +587,6 @@ template <class ScalarT, class TraitsT>
         return detail::intersects_any(r, caps[0], caps[1], boundary_no_caps_of(c));
     }
 }
-
-// ray - tube
 template <class ScalarT>
 [[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, cylinder_boundary_no_caps<3, ScalarT> const& c)
 {
@@ -638,6 +636,27 @@ template <class ScalarT>
     return {hits, hit_cnt};
 }
 
+// ray - inf_cylinder
+template <class ScalarT>
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<2, ScalarT> const& r, inf_cylinder_boundary<2, ScalarT> const& c)
+{
+    const auto n = perpendicular(c.axis.dir);
+    const auto d = dot(r.dir, n);
+    if (d == ScalarT(0)) // ray parallel to inf_cylinder
+        return {};
+
+    const auto dist = dot(c.axis.pos - r.origin, n);
+    const auto [tMin, tMax] = minmax((dist - c.radius) / d, (dist + c.radius) / d);
+
+    if (tMin >= ScalarT(0))
+        return {tMin, tMax};
+
+    if (tMax >= ScalarT(0))
+        return tMax;
+
+    return {};
+}
+
 // ray - triangle2
 template <class ScalarT>
 [[nodiscard]] constexpr optional<ray_interval<ScalarT>> intersection_parameter(ray<2, ScalarT> const& r, triangle<2, ScalarT> const& t)
@@ -725,8 +744,7 @@ template <class ScalarT>
     const auto t1 = (-b - sqrtD) / (ScalarT(2) * a);
     const auto t2 = (-b + sqrtD) / (ScalarT(2) * a);
 
-    auto tMin = min(t1, t2);
-    auto tMax = max(t1, t2);
+    const auto [tMin, tMax] = minmax(t1, t2);
 
     if (tMin >= ScalarT(0))
         return {tMin, tMax};
diff --git a/src/typed-geometry/types/objects/inf_cylinder.hh b/src/typed-geometry/types/objects/inf_cylinder.hh
index 2a97962..f4c1a16 100644
--- a/src/typed-geometry/types/objects/inf_cylinder.hh
+++ b/src/typed-geometry/types/objects/inf_cylinder.hh
@@ -29,6 +29,7 @@ using inf_cylinder_boundary = inf_cylinder<D, ScalarT, boundary_tag>;
 template <int D, class ScalarT, class TraitsT>
 struct inf_cylinder
 {
+    using scalar_t = ScalarT;
     using pos_t = pos<D, ScalarT>;
     using dir_t = dir<D, ScalarT>;
     using line_t = line<D, ScalarT>;
-- 
GitLab


From 49b66510f1807a1457dce468f8e63a4ba6a451ee Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Tue, 17 Nov 2020 22:12:06 +0100
Subject: [PATCH 29/42] Added ::scalar_t to line and ray

---
 src/typed-geometry/types/objects/line.hh | 6 +++---
 src/typed-geometry/types/objects/ray.hh  | 6 +++---
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/typed-geometry/types/objects/line.hh b/src/typed-geometry/types/objects/line.hh
index d0f84b5..02d2891 100644
--- a/src/typed-geometry/types/objects/line.hh
+++ b/src/typed-geometry/types/objects/line.hh
@@ -3,7 +3,6 @@
 #include <typed-geometry/types/scalars/default.hh>
 #include "../dir.hh"
 #include "../pos.hh"
-#include "../size.hh"
 #include "traits.hh"
 
 // A line has a direction and a point lying on it
@@ -46,8 +45,9 @@ using uline4 = line<4, u32>;
 template <int D, class ScalarT>
 struct line
 {
-    using dir_t = tg::dir<D, ScalarT>;
-    using pos_t = tg::pos<D, ScalarT>;
+    using scalar_t = ScalarT;
+    using dir_t = dir<D, ScalarT>;
+    using pos_t = pos<D, ScalarT>;
 
     pos_t pos;
     dir_t dir;
diff --git a/src/typed-geometry/types/objects/ray.hh b/src/typed-geometry/types/objects/ray.hh
index 402f3d0..a376315 100644
--- a/src/typed-geometry/types/objects/ray.hh
+++ b/src/typed-geometry/types/objects/ray.hh
@@ -3,7 +3,6 @@
 #include <typed-geometry/types/scalars/default.hh>
 #include "../dir.hh"
 #include "../pos.hh"
-#include "../size.hh"
 #include "traits.hh"
 
 namespace tg
@@ -44,8 +43,9 @@ using uray4 = ray<4, u32>;
 template <int D, class ScalarT>
 struct ray
 {
-    using dir_t = tg::dir<D, ScalarT>;
-    using pos_t = tg::pos<D, ScalarT>;
+    using scalar_t = ScalarT;
+    using dir_t = dir<D, ScalarT>;
+    using pos_t = pos<D, ScalarT>;
 
     pos_t origin;
     dir_t dir;
-- 
GitLab


From 1b11ecffb58cda9dd6ac80da1dae8da63e104809 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Wed, 18 Nov 2020 22:30:24 +0100
Subject: [PATCH 30/42] Changed all explicit intersections from ray to line
 intersections and added default conversions from line to ray

---
 src/typed-geometry/detail/utility.hh          |  12 +
 .../functions/objects/intersection.hh         | 550 +++++++++---------
 2 files changed, 292 insertions(+), 270 deletions(-)

diff --git a/src/typed-geometry/detail/utility.hh b/src/typed-geometry/detail/utility.hh
index 497813f..2f0f579 100644
--- a/src/typed-geometry/detail/utility.hh
+++ b/src/typed-geometry/detail/utility.hh
@@ -35,7 +35,19 @@ template <>
 struct priority_tag<0>
 {
 };
+
+// see https://stackoverflow.com/questions/44395169/why-is-sfinae-on-if-constexpr-not-allowed
+template <template <class...> class, class, class...>
+struct can_apply : std::false_type
+{
+};
+template <template <class...> class Z, class... Ts>
+struct can_apply<Z, std::void_t<Z<Ts...>>, Ts...> : std::true_type
+{
+};
 } // namespace detail
+template <template <class...> class Z, class... Ts>
+constexpr bool can_apply = detail::can_apply<Z, void, Ts...>::value;
 
 template <class...>
 constexpr bool always_false = false;
diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index dd25a67..b25bcc4 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -32,7 +32,7 @@
 // intersects(a, b)              -> bool
 // intersection(a, b)            -> ???
 // intersection_safe(a, b)       -> optional<???>
-// intersection_parameter(a, b)  -> coords? (for a ray: ray_hits<N, ScalarT> or optional<ray_interval> (when b is solid))
+// intersection_parameter(a, b)  -> coords? (for a line or a ray: ray_hits<N, ScalarT> or optional<ray_interval> (when b is solid))
 // intersection_parameters(a, b) -> pair<coords, coords>?
 // intersection_exact(a, b)      -> variant
 // closest_intersection(a, b)            -> position/object (for a ray: optional<pos>)
@@ -52,15 +52,22 @@
 //  - for more elaborate ray-tracing, a future ray_cast function will exist (which also returns the intersection normal)
 
 // Implementation guidelines:
-// explicit intersection_parameter(ray, obj), which gives closest_intersection_parameter and intersection position(s)
-//          intersection_parameter(ray, obj_boundary) already gives intersection_parameter(ray, obj) when boundary_of(obj) is defined and object is convex and finite
-//          for infinite objects, a separate solid version is needed when the boundary version does not distinguish between completely contained and completely outside
-// explicit closest_intersection_parameter(ray, obj) if this is faster than computing all intersections
+// 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
+//
 // explicit intersects(obj, aabb), which gives intersects(aabb, obj)
 //
 // for convex compound objects (like cylinder or pyramids), decompose the object into primitive shapes and pass them to a helper function:
-// - call merge_hits(ray, objPart1, objPart2, ...) in the implementation of intersection_parameter(ray, obj_boundary)
-// - call intersects_any(ray, objPart1, objPart2, ...) in the implementation of intersects(ray, obj<TraitsT>), which can shortcut and be faster than the default
+// - 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
 
 
 namespace tg
@@ -77,6 +84,7 @@ struct ray_hits
 {
     static constexpr bool is_ray_hits = true; // tag
     static constexpr int max_hits = MaxHits;
+    using hit_t = HitT;
 
     template <class OtherT>
     using as_ray_hits = ray_hits<MaxHits, OtherT>;
@@ -94,6 +102,11 @@ struct ray_hits
         TG_ASSERT(_size > 0);
         return _hit[0];
     }
+    [[nodiscard]] HitT const& last() const
+    {
+        TG_ASSERT(_size > 0);
+        return _hit[_size - 1];
+    }
 
     [[nodiscard]] HitT const* begin() const { return _hit; }
     [[nodiscard]] HitT const* end() const { return _hit + _size; }
@@ -118,28 +131,72 @@ private:
 template <class ScalarT>
 struct ray_interval
 {
+    static constexpr bool is_ray_interval = true; // tag
+
     ScalarT start;
     ScalarT end;
 
-    [[nodiscard]] constexpr bool is_unbounded() const { return end == tg::max<ScalarT>(); }
+    [[nodiscard]] constexpr bool is_unbounded() const { return end == tg::max<ScalarT>() || start == tg::min<ScalarT>(); }
 };
 
 
+// ====================================== 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>
+[[nodiscard]] constexpr ray_hits<2, ScalarT> merge_hits(line<D, ScalarT> const& line, Objs const&... objs)
+{
+    ScalarT hits[2];
+    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};
+}
+
+// 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) || ...);
+}
+
+template <class A, class B>
+using try_closest_intersection_parameter = decltype(closest_intersection_parameter(std::declval<A const&>(), std::declval<B const&>()).has_value());
+}
+
+
 // ====================================== Default Implementations ======================================
 // TODO: intersection_parameter from intersection_parameters
 
 // returns whether two objects intersect
+// if intersection is available and applicable, use that
 template <class A, class B>
-[[nodiscard]] constexpr auto intersects(A const& a, B const& b) -> decltype(intersection(a, b).has_value())
+[[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())>
 {
     return intersection(a, b).has_value();
 }
 
 // if closest intersection parameter is available and applicable, use that
-template <int D, class ScalarT, class Obj>
-[[nodiscard]] constexpr auto intersects(ray<D, ScalarT> const& r, Obj const& obj) -> decltype(closest_intersection_parameter(r, obj).has_value())
+template <class A, class B>
+[[nodiscard]] constexpr auto intersects(A const& a, B const& b) -> decltype(closest_intersection_parameter(a, b).has_value())
 {
-    return closest_intersection_parameter(r, obj).has_value();
+    return closest_intersection_parameter(a, b).has_value();
 }
 
 // parameters for intersects with aabb can switch order
@@ -188,45 +245,44 @@ template <class A, class B>
 }
 
 // if ray_hits intersection parameter is available, use that
-template <int D, class ScalarT, class Obj>
-[[nodiscard]] constexpr auto intersection(ray<D, ScalarT> const& r, Obj const& obj) ->
-    typename decltype(intersection_parameter(r, obj))::template as_ray_hits<pos<D, ScalarT>>
+template <class A, class B>
+[[nodiscard]] constexpr auto intersection(A const& a, B const& b) ->
+    typename decltype(intersection_parameter(a, b))::template as_ray_hits<typename A::pos_t>
 {
-    auto ts = intersection_parameter(r, obj);
-    pos<D, ScalarT> hits[ts.max_hits];
+    auto ts = intersection_parameter(a, b);
+    typename A::pos_t hits[ts.max_hits];
     for (auto i = 0; i < ts.size(); ++i)
-        hits[i] = r[ts[i]];
+        hits[i] = a[ts[i]];
     return {hits, ts.size()};
 }
 
 // if an optional ray_interval intersection parameter is available, use that
-template <int D, class ScalarT, class Obj, typename = std::enable_if_t<std::is_same_v<decltype(
-    intersection_parameter(std::declval<ray<D, ScalarT>>(), std::declval<Obj>())), optional<ray_interval<ScalarT>>>>>
-[[nodiscard]] constexpr optional<segment<D, ScalarT>> intersection(ray<D, ScalarT> const& r, Obj const& obj)
+template <class A, class B>
+[[nodiscard]] constexpr auto intersection(A const& a, B const& b)
+    -> std::enable_if<decltype(intersection_parameter(a, b).value())::is_ray_interval, // condition
+    optional<segment<object_traits<A>::domain_dimension, typename object_traits<A>::scalar_t>>> // return type
 {
-    auto ts = intersection_parameter(r, obj);
+    auto ts = intersection_parameter(a, b);
     if (ts.has_value())
-        return {{r[ts.value().start], r[ts.value().end]}};
+        return {{a[ts.value().start], a[ts.value().end]}};
     return {};
 }
 
 // if ray_hits intersection parameter is available, use that
-template <int D, class ScalarT, class Obj>
-[[nodiscard]] constexpr auto closest_intersection_parameter(ray<D, ScalarT> const& r, Obj const& obj)
-    -> enable_if<decltype(intersection_parameter(r, obj))::is_ray_hits, optional<ScalarT>>
+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>
 {
-    auto hits = intersection_parameter(r, obj);
+    const auto hits = intersection_parameter(a, b);
     if (hits.any())
         return hits.first();
     return {};
 }
 
 // if ray_interval intersection parameter is available, use that
-template <int D, class ScalarT, class Obj>
-[[nodiscard]] constexpr auto closest_intersection_parameter(ray<D, ScalarT> const& r, Obj const& obj)
-    -> enable_if<std::is_same_v<decltype(intersection_parameter(r, obj)), optional<ray_interval<ScalarT>>>, optional<ScalarT>>
+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)>
 {
-    auto hits = intersection_parameter(r, obj);
+    const auto hits = intersection_parameter(a, b);
     if (hits.has_value())
         return hits.value().start;
     return {};
@@ -234,22 +290,33 @@ template <int D, class ScalarT, class Obj>
 
 // if boundary_of a solid object returns ray_hits, use this to construct the ray_interval result of the solid intersection
 template <int D, class ScalarT, class Obj>
-[[nodiscard]] constexpr auto intersection_parameter(ray<D, ScalarT> const& r, Obj const& obj)
+[[nodiscard]] constexpr auto intersection_parameter(line<D, ScalarT> const& l, Obj const& obj)
     -> enable_if<!std::is_same_v<Obj, decltype(boundary_of(obj))>, optional<ray_interval<ScalarT>>>
 {
-    const ray_hits<2, ScalarT> inter = intersection_parameter(r, boundary_of(obj));
+    const ray_hits<2, ScalarT> inter = intersection_parameter(l, boundary_of(obj));
 
-    if (inter.size() == 0)
-        return {};
+    if (inter.size() == 2)
+        return {{inter[0], inter[1]}};
 
-    if (inter.size() == 1)
+    if constexpr (object_traits<Obj>::is_finite)
     {
-        TG_ASSERT(contains(obj, r.origin));
-        return {{typename Obj::scalar_t(0), inter.first()}};
+        TG_ASSERT(inter.size() == 0); // a line intersects a finite solid object either twice or not at all
+        return {};
     }
+    else
+    {
+        if (inter.size() == 0)
+        {
+            if (contains(obj, l.pos))
+                return {{tg::min<ScalarT>(), tg::max<ScalarT>()}};
+            return {};
+        }
 
-    TG_ASSERT(inter.size() == 2);
-    return {{inter[0], inter[1]}};
+        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()}};
+    }
 }
 
 // intersection between point and obj is same as contains
@@ -271,131 +338,117 @@ constexpr optional<pos<D, ScalarT>> intersection(Obj const& obj, pos<D, ScalarT>
 }
 
 
-// ====================================== Helper functions ======================================
+// ====================================== Ray Intersections from Line Intersections ======================================
 
-namespace detail
-{
-// intersects the given ray with all given objects and returns the concatenated intersections. A maximal number of 2 intersections is assumed.
-template <int D, class ScalarT, class... ObjsT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> merge_hits(ray<D, ScalarT> const& ray, ObjsT const&... objs)
+template <int D, class ScalarT, class Obj>
+[[nodiscard]] constexpr auto intersection_parameter(ray<D, ScalarT> const& ray, Obj const& obj)
+    -> decltype(intersection_parameter(inf_of(ray), obj).is_ray_hits, intersection_parameter(inf_of(ray), obj))
 {
-    ScalarT hits[2];
-    hits[0] = tg::max<ScalarT>();
-    hits[1] = tg::min<ScalarT>();
-    auto numHits = 0;
+    const auto inter = intersection_parameter(inf_of(ray), obj);
+    constexpr auto maxHits = inter.max_hits;
+    if (!inter.any() || inter.last() < ScalarT(0))
+        return {};
 
-    const auto find_hits = [&](const auto& obj) {
-        const auto inters = intersection_parameter(ray, obj);
-        for (const auto& inter : inters)
-        {
-            hits[0] = tg::min(hits[0], inter);
-            hits[1] = tg::max(hits[1], inter);
-            numHits++;
-        }
-    };
-    (find_hits(objs), ...);
+    if constexpr (maxHits == 2)
+    {
+        if (inter.first() < ScalarT(0))
+            return {inter[1]};
+    }
+    else
+        static_assert(maxHits == 1, "Only up to two intersections implemented");
 
-    TG_ASSERT(numHits <= 2);
-    return {hits, numHits};
+    return inter;
 }
-
-// returns true, iff the given ray intersects any of the given objects (with short-circuiting after the first intersection)
-template <int D, class ScalarT, class... ObjsT>
-[[nodiscard]] constexpr bool intersects_any(ray<D, ScalarT> const& ray, ObjsT const&... objs)
+template <int D, class ScalarT, class Obj>
+[[nodiscard]] constexpr auto intersection_parameter(ray<D, ScalarT> const& ray, Obj const& obj)
+    -> decltype(intersection_parameter(inf_of(ray), obj).value().is_ray_interval, intersection_parameter(inf_of(ray), obj))
 {
-    return (intersects(ray, objs) || ...);
-}
+    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;
 }
 
-// ====================================== Ray - Object Intersections ======================================
+// ====================================== Line - Object Intersections ======================================
 
-// ray - point
+// line - point
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<1, ScalarT> const& r, pos<1, ScalarT> const& p)
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<1, ScalarT> const& l, pos<1, ScalarT> const& p)
 {
-    const auto t = coordinates(r, p);
-    if (t >= ScalarT(0))
-        return t;
-    return {};
+    return coordinates(l, p);
 }
 
-// ray - line
+// line - line
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<2, ScalarT> const& r, line<2, ScalarT> const& l)
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<2, ScalarT> const& l0, line<2, ScalarT> const& l1)
 {
-    // r.origin + r.dir * t.x == l.pos + l.dir * t.y  <=>  (r.dir | -l.dir) * (t.x | t.y)^T == l.pos - r.origin
-    auto M = mat<2, 2, ScalarT>::from_cols(r.dir, -l.dir);
-    auto t = inverse(M) * (l.pos - r.origin);
-    if (t.x < ScalarT(0))
-        return {};
+    // 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);
     return t.x;
 }
 
-// ray - ray
+// line - ray
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<2, ScalarT> const& r0, ray<2, ScalarT> const& r1)
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<2, ScalarT> const& l, ray<2, ScalarT> const& r)
 {
-    auto M = mat<2, 2, ScalarT>::from_cols(r0.dir, -r1.dir);
-    auto t = inverse(M) * (r1.origin - r0.origin);
-    if (t.x < ScalarT(0) || t.y < ScalarT(0))
+    auto M = mat<2, 2, ScalarT>::from_cols(l.dir, -r.dir);
+    auto t = inverse(M) * (r.origin - l.pos);
+    if (t.y < ScalarT(0))
         return {};
     return t.x;
 }
 
-// ray - segment
+// line - segment
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<2, ScalarT> const& r, segment<2, ScalarT> const& s)
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<2, ScalarT> const& l, segment<2, ScalarT> const& s)
 {
-    auto M = mat<2, 2, ScalarT>::from_cols(r.dir, s.pos0 - s.pos1);
-    auto t = inverse(M) * (s.pos0 - r.origin);
-    if (t.x < ScalarT(0) || t.y < ScalarT(0) || t.y > ScalarT(1))
+    auto M = mat<2, 2, ScalarT>::from_cols(l.dir, s.pos0 - s.pos1);
+    auto t = inverse(M) * (s.pos0 - l.pos);
+    if (t.y < ScalarT(0) || t.y > ScalarT(1))
         return {};
     return t.x;
 }
 
-// ray - plane
+// line - plane
 template <int D, class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, plane<D, ScalarT> const& p)
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<D, ScalarT> const& l, plane<D, ScalarT> const& p)
 {
-    // if plane normal and ray direction are orthogonal there is no intersection
-    const auto dotND = dot(p.normal, r.dir);
-    if (dotND == ScalarT(0))
+    const auto dotND = dot(p.normal, l.dir);
+    if (dotND == ScalarT(0)) // if plane normal and line direction are orthogonal, there is no intersection
         return {};
 
-    // plane: <x, p.normal> = p.dis
-    // ray: x = r.origin + t * r.dir
-    // =>
-    // <r.origin + t * r.dir, p.normal> = p.dis
-    // t = (p.dis - <r.origin, p.normal>) / <r.dir, p.normal>
-
-    const auto t = (p.dis - dot(p.normal, r.origin)) / dotND;
-
-    // check whether plane lies behind ray
-    if (t < ScalarT(0))
-        return {};
-
-    return t;
+    // <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;
 }
 
-// ray - halfspace
+// line - halfspace
 template <int D, class ScalarT>
-[[nodiscard]] constexpr optional<ray_interval<ScalarT>> intersection_parameter(ray<D, ScalarT> const& r, halfspace<D, ScalarT> const& h)
+[[nodiscard]] constexpr optional<ray_interval<ScalarT>> intersection_parameter(line<D, ScalarT> const& l, halfspace<D, ScalarT> const& h)
 {
-    // check if ray origin is already contained in the halfspace
-    const auto dotND = dot(h.normal, r.dir);
-    const auto dist = signed_distance(r.origin, h);
-    if (dist <= ScalarT(0))
+    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
     {
-        if (dotND <= ScalarT(0)) // if ray point away from the surface, the entire ray intersects
-            return {{ScalarT(0), tg::max<ScalarT>()}};
-        return {{ScalarT(0), -dist / dotND}}; // otherwise the segment from ray origin until surface intersection is contained
+        if (dist <= ScalarT(0))
+            return {{tg::min<ScalarT>(), tg::max<ScalarT>()}}; // completely contained
+        return {}; // completely outside
     }
 
-    // if ray points away from the halfspace there is no intersection
-    if (dotND >= ScalarT(0))
-        return {};
-
-    return {{-dist / dotND, tg::max<ScalarT>()}};
+    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
 }
 template <int D, class ScalarT>
 [[nodiscard]] constexpr optional<ScalarT> closest_intersection_parameter(ray<D, ScalarT> const& r, halfspace<D, ScalarT> const& h)
@@ -413,46 +466,43 @@ template <int D, class ScalarT>
     return -dist / dotND;
 }
 
-// ray - aabb
+// line - aabb
 template <int D, class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, aabb_boundary<D, ScalarT> const& b)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, aabb_boundary<D, ScalarT> const& b)
 {
     auto tFirst = tg::min<ScalarT>();
     auto tSecond = tg::max<ScalarT>();
     for (auto i = 0; i < D; ++i)
     {
-        if (abs(r.dir[i]) > ScalarT(100) * tg::epsilon<ScalarT>)
+        if (abs(l.dir[i]) > ScalarT(100) * tg::epsilon<ScalarT>)
         {
-            const auto tMin = (b.min[i] - r.origin[i]) / r.dir[i];
-            const auto tMax = (b.max[i] - r.origin[i]) / r.dir[i];
+            const auto tMin = (b.min[i] - l.pos[i]) / l.dir[i];
+            const auto tMax = (b.max[i] - l.pos[i]) / l.dir[i];
             tFirst = max(tFirst, min(tMin, tMax));
             tSecond = min(tSecond, max(tMin, tMax));
         }
-        else if (r.origin[i] < b.min[i] || r.origin[i] > b.max[i])
+        else if (l.pos[i] < b.min[i] || l.pos[i] > b.max[i])
             return {}; // ray parallel to this axis and outside of the aabb
     }
 
-    // no intersection if ray starts behind the aabb or if it misses the aabb
-    if (tSecond < ScalarT(0) || tFirst > tSecond)
+    // no intersection if line misses the aabb
+    if (tFirst > tSecond)
         return {};
 
-    if (tFirst < ScalarT(0))
-        return tSecond;
-
     return {tFirst, tSecond};
 }
 
-// ray - box
+// line - box
 template <int D, class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, box_boundary<D, ScalarT> const& b)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, box_boundary<D, ScalarT> const& b)
 {
-    const auto bMin = b[comp<D, ScalarT>(-1)] - r.origin;
-    const auto bMax = b[comp<D, ScalarT>(1)] - r.origin;
+    const auto bMin = b[comp<D, ScalarT>(-1)] - l.pos;
+    const auto bMax = b[comp<D, ScalarT>(1)] - l.pos;
     auto tFirst = tg::min<ScalarT>();
     auto tSecond = tg::max<ScalarT>();
     for (auto i = 0; i < D; ++i)
     {
-        const auto rDir = dot(r.dir, b.half_extents[i]);
+        const auto rDir = dot(l.dir, b.half_extents[i]);
         if (abs(rDir) > ScalarT(100) * tg::epsilon<ScalarT>)
         {
             const auto tMin = dot(bMin, b.half_extents[i]) / rDir;
@@ -464,161 +514,154 @@ template <int D, class ScalarT>
             return {}; // ray parallel to this half_extents axis and outside of the box
     }
 
-    // no intersection if ray starts behind the box or if it misses the box
-    if (tSecond < ScalarT(0) || tFirst > tSecond)
+    // no intersection if line misses the box
+    if (tFirst > tSecond)
         return {};
 
-    if (tFirst < ScalarT(0))
-        return tSecond;
-
     return {tFirst, tSecond};
 }
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, box<2, ScalarT, 3> const& b)
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<3, ScalarT> const& l, box<2, ScalarT, 3> const& b)
 {
-    const auto t = intersection_parameter(r, plane<3, ScalarT>(normal_of(b), b.center));
-    if (!t.any())
+    const auto t = intersection_parameter(l, plane<3, ScalarT>(normal_of(b), b.center));
+    if (!t.any()) // line parallel to box plane
         return {};
 
-    const auto p = r[t.first()] - b.center;
+    const auto p = l[t.first()] - b.center;
     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;
 }
 
-// ray - disk
+// line - disk
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, sphere<2, ScalarT, 3> const& d)
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<3, ScalarT> const& l, sphere<2, ScalarT, 3> const& d)
 {
-    const auto t = intersection_parameter(r, plane<3, ScalarT>(d.normal, d.center));
-    if (!t.any())
+    const auto t = intersection_parameter(l, plane<3, ScalarT>(d.normal, d.center));
+    if (!t.any()) // line parallel to disk plane
         return {};
 
-    const auto p = r[t.first()];
+    const auto p = l[t.first()];
     if (distance_sqr(p, d.center) > d.radius * d.radius)
         return {};
 
     return t;
 }
 
-// ray - sphere
+// line - sphere
 template <int D, class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, sphere_boundary<D, ScalarT> const& s)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, sphere_boundary<D, ScalarT> const& s)
 {
-    auto t = dot(s.center - r.origin, r.dir);
+    auto t = dot(s.center - l.pos, l.dir);
 
-    auto d_sqr = distance_sqr(r[t], s.center);
+    auto d_sqr = distance_sqr(l[t], s.center);
     auto r_sqr = s.radius * s.radius;
     if (d_sqr > r_sqr)
         return {};
 
     auto dt = sqrt(r_sqr - d_sqr);
-
-    if (t - dt >= 0)
-        return {t - dt, t + dt};
-
-    if (t + dt >= 0)
-        return {t + dt};
-
-    return {};
+    return {t - dt, t + dt};
 }
 
-// ray - hemisphere
+// line - hemisphere
 template <int D, class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, hemisphere_boundary<D, ScalarT> const& h)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, hemisphere_boundary<D, ScalarT> const& h)
 {
-    return detail::merge_hits(r, caps_of(h), boundary_no_caps_of(h));
+    return detail::merge_hits(l, caps_of(h), boundary_no_caps_of(h));
 }
 template <int D, class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<D, ScalarT> const& r, hemisphere_boundary_no_caps<D, ScalarT> const& h)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, hemisphere_boundary_no_caps<D, ScalarT> const& h)
 {
     ScalarT hits[2];
     auto numHits = 0;
-    const auto sphereHits = intersection_parameter(r, sphere_boundary<D, ScalarT>(h.center, h.radius));
+    const auto sphereHits = intersection_parameter(l, sphere_boundary<D, ScalarT>(h.center, h.radius));
     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)
-        if (contains(halfSpace, r[hit]))
+        if (contains(halfSpace, l[hit]))
             hits[numHits++] = hit;
 
     return {hits, numHits};
 }
-template <int D, class ScalarT, class TraitsT>
-[[nodiscard]] constexpr bool intersects(ray<D, ScalarT> const& r, hemisphere<D, ScalarT, TraitsT> const& h)
+template <class Obj, int D, class ScalarT, class TraitsT>
+[[nodiscard]] constexpr bool intersects(Obj const& lr, hemisphere<D, ScalarT, TraitsT> const& h)
 {
+    static_assert(object_traits<Obj>::is_infinite, "For finite objects, complete containment within boundary has to be considered as well");
     if constexpr (std::is_same_v<TraitsT, boundary_no_caps_tag>)
-        return intersection_parameter(r, h).any();
+        return intersection_parameter(lr, h).any();
     else
-        return detail::intersects_any(r, caps_of(h), boundary_no_caps_of(h));
+        return detail::intersects_any(lr, caps_of(h), boundary_no_caps_of(h));
 }
 
-// ray - capsule
+// line - capsule
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, capsule_boundary<3, ScalarT> const& c)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, capsule_boundary<3, ScalarT> const& c)
 {
     using caps_t = hemisphere_boundary_no_caps<3, ScalarT>;
     const auto n = direction(c);
-    return detail::merge_hits(r, caps_t(c.axis.pos0, c.radius, -n), caps_t(c.axis.pos1, c.radius, n),
+    return detail::merge_hits(l, caps_t(c.axis.pos0, c.radius, -n), caps_t(c.axis.pos1, c.radius, n),
                               cylinder_boundary_no_caps<3, ScalarT>(c.axis, c.radius));
 }
-template <class ScalarT, class TraitsT>
-[[nodiscard]] constexpr bool intersects(ray<3, ScalarT> const& r, capsule<3, ScalarT, TraitsT> const& c)
+template <class Obj, class ScalarT, class TraitsT>
+[[nodiscard]] constexpr bool intersects(Obj const& lr, capsule<3, ScalarT, TraitsT> const& c)
 {
+    static_assert(object_traits<Obj>::is_infinite, "For finite objects, complete containment within boundary has to be considered as well");
     using caps_t = sphere_boundary<3, ScalarT>; // spheres are faster than hemispheres and equivalent for the yes/no decision
-    return detail::intersects_any(r, caps_t(c.axis.pos0, c.radius), caps_t(c.axis.pos1, c.radius),
+    return detail::intersects_any(lr, caps_t(c.axis.pos0, c.radius), caps_t(c.axis.pos1, c.radius),
                                   cylinder_boundary_no_caps<3, ScalarT>(c.axis, c.radius));
 }
 
-// ray - cylinder
+// line - cylinder
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, cylinder_boundary<3, ScalarT> const& c)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, cylinder_boundary<3, ScalarT> const& c)
 {
     const auto caps = caps_of(c);
-    return detail::merge_hits(r, caps[0], caps[1], boundary_no_caps_of(c));
+    return detail::merge_hits(l, caps[0], caps[1], boundary_no_caps_of(c));
 }
-template <class ScalarT, class TraitsT>
-[[nodiscard]] constexpr bool intersects(ray<3, ScalarT> const& r, cylinder<3, ScalarT, TraitsT> const& c)
+template <class Obj, class ScalarT, class TraitsT>
+[[nodiscard]] constexpr bool intersects(Obj const& lr, cylinder<3, ScalarT, TraitsT> const& c)
 {
+    static_assert(object_traits<Obj>::is_infinite, "For finite objects, complete containment within boundary has to be considered as well");
     if constexpr (std::is_same_v<TraitsT, boundary_no_caps_tag>)
-        return intersection_parameter(r, c).any();
+        return intersection_parameter(lr, c).any();
     else
     {
         const auto caps = caps_of(c);
-        return detail::intersects_any(r, caps[0], caps[1], boundary_no_caps_of(c));
+        return detail::intersects_any(lr, caps[0], caps[1], boundary_no_caps_of(c));
     }
 }
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, cylinder_boundary_no_caps<3, ScalarT> const& c)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, cylinder_boundary_no_caps<3, ScalarT> const& c)
 {
     auto cdir = direction(c);
-    auto cosA = dot(cdir, r.dir);
+    auto cosA = dot(cdir, l.dir);
     auto sinA_sqr = 1 - cosA * cosA;
 
     if (sinA_sqr <= 0)
-        return {};
-
-    // compute closest points of the two lines
-    auto origDiff = r.origin - c.axis.pos0;
-    auto fRay = dot(r.dir, origDiff);
-    auto fLine = dot(cdir, origDiff);
-    auto tRay = (cosA * fLine - fRay) / sinA_sqr;
-    auto tLine = (fLine - cosA * fRay) / sinA_sqr;
-
-    auto closest_on_ray = r.origin + tRay * r.dir;
-    auto closest_on_line = c.axis.pos0 + tLine * cdir;
-    auto line_ray_dist_sqr = distance_sqr(closest_on_ray, closest_on_line);
+        return {}; // line and cylinder are parallel
+
+    // compute closest points of line and cylinder axis
+    auto origDiff = l.pos - c.axis.pos0;
+    auto fLine = dot(l.dir, origDiff);
+    auto fAxis = dot(cdir, origDiff);
+    auto tLine = (cosA * fAxis - fLine) / sinA_sqr;
+    auto tAxis = (fAxis - cosA * fLine) / sinA_sqr;
+
+    auto closest_on_line = l.pos + tLine * l.dir;
+    auto closest_on_axis = c.axis.pos0 + tAxis * cdir;
+    auto line_axis_dist_sqr = distance_sqr(closest_on_line, closest_on_axis);
     auto cyl_radius_sqr = c.radius * c.radius;
 
-    if (line_ray_dist_sqr > cyl_radius_sqr)
-        return {};
+    if (line_axis_dist_sqr > cyl_radius_sqr)
+        return {}; // line misses the cylinder
 
     // radius in 2D slice
-    auto r_2D = sqrt(cyl_radius_sqr - line_ray_dist_sqr);
+    auto r_2D = sqrt(cyl_radius_sqr - line_axis_dist_sqr);
 
     // infinite tube intersection
     auto s = r_2D / sqrt(sinA_sqr);
-    auto cyl_intersection_0 = closest_on_ray - s * r.dir;
-    auto cyl_intersection_1 = closest_on_ray + s * r.dir;
+    auto cyl_intersection_0 = closest_on_line - s * l.dir;
+    auto cyl_intersection_1 = closest_on_line + s * l.dir;
 
     // project onto line segment
     auto line_length = length(c.axis);
@@ -628,45 +671,38 @@ template <class ScalarT>
     ScalarT hits[2];
     int hit_cnt = 0;
 
-    if (tRay - s >= 0 && 0 <= lambda_0 && lambda_0 < line_length)
-        hits[hit_cnt++] = tRay - s;
-    if (tRay + s >= 0 && 0 <= lambda_1 && lambda_1 < line_length)
-        hits[hit_cnt++] = tRay + s;
+    if (0 <= lambda_0 && lambda_0 < line_length)
+        hits[hit_cnt++] = tLine - s;
+    if (0 <= lambda_1 && lambda_1 < line_length)
+        hits[hit_cnt++] = tLine + s;
 
     return {hits, hit_cnt};
 }
 
-// ray - inf_cylinder
+// line - inf_cylinder
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<2, ScalarT> const& r, inf_cylinder_boundary<2, ScalarT> const& c)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<2, ScalarT> const& l, inf_cylinder_boundary<2, ScalarT> const& c)
 {
     const auto n = perpendicular(c.axis.dir);
-    const auto d = dot(r.dir, n);
-    if (d == ScalarT(0)) // ray parallel to inf_cylinder
+    const auto d = dot(l.dir, n);
+    if (d == ScalarT(0)) // line parallel to inf_cylinder
         return {};
 
-    const auto dist = dot(c.axis.pos - r.origin, n);
+    const auto dist = dot(c.axis.pos - l.pos, n);
     const auto [tMin, tMax] = minmax((dist - c.radius) / d, (dist + c.radius) / d);
-
-    if (tMin >= ScalarT(0))
-        return {tMin, tMax};
-
-    if (tMax >= ScalarT(0))
-        return tMax;
-
-    return {};
+    return {tMin, tMax};
 }
 
-// ray - triangle2
+// line - triangle2
 template <class ScalarT>
-[[nodiscard]] constexpr optional<ray_interval<ScalarT>> intersection_parameter(ray<2, ScalarT> const& r, triangle<2, ScalarT> const& t)
+[[nodiscard]] constexpr optional<ray_interval<ScalarT>> intersection_parameter(line<2, ScalarT> const& l, triangle<2, ScalarT> const& t)
 {
     ScalarT closestIntersection = tg::max<ScalarT>();
     ScalarT furtherIntersection = tg::min<ScalarT>();
     auto numIntersections = 0;
     for (const auto& edge : edges_of(t))
     {
-        const auto inter = intersection_parameter(r, edge);
+        const auto inter = intersection_parameter(l, edge);
         if (inter.any())
         {
             numIntersections++;
@@ -677,62 +713,57 @@ template <class ScalarT>
 
     if (numIntersections == 0)
         return {}; // No intersection
-    if (numIntersections == 1)
-    {
-        // ray started within the triangle
-        TG_ASSERT(contains(t, r.origin));
-        return {{ScalarT(0), furtherIntersection}};
-    }
+
+    TG_ASSERT(numIntersections == 2);
     return {{closestIntersection, furtherIntersection}};
 }
 
-// ray - triangle3
+// line - triangle3
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(ray<3, ScalarT> const& r,
+[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<3, ScalarT> const& l,
                                                                  triangle<3, ScalarT> const& t,
                                                                  dont_deduce<ScalarT> eps = 100 * tg::epsilon<ScalarT>)
 {
     auto e1 = t.pos1 - t.pos0;
     auto e2 = t.pos2 - t.pos0;
 
-    auto pvec = tg::cross(r.dir, e2);
+    auto pvec = tg::cross(l.dir, e2);
     auto det = dot(pvec, e1);
 
     if (det < ScalarT(0))
     {
         std::swap(e1, e2);
-        pvec = tg::cross(r.dir, e2);
+        pvec = tg::cross(l.dir, e2);
         det = -det;
     }
 
     if (det < eps)
         return {};
 
-    auto tvec = r.origin - t.pos0;
+    auto tvec = l.pos - t.pos0;
     auto u = dot(tvec, pvec);
     if (u < ScalarT(0) || u > det)
         return {};
 
     auto qvec = cross(tvec, e1);
-    auto v = dot(r.dir, qvec);
+    auto v = dot(l.dir, qvec);
     if (v < ScalarT(0) || v + u > det)
         return {};
 
-    auto lambda = (ScalarT(1) / det) * dot(e2, qvec);
-    return (lambda > ScalarT(0)) ? lambda : ray_hits<1, ScalarT>();
+    return (ScalarT(1) / det) * dot(e2, qvec);
 }
 
-// ray - quadric_boundary (as an isosurface, not error quadric)
+// line - quadric_boundary (as an isosurface, not error quadric)
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, quadric<3, ScalarT> const& q)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, quadric<3, ScalarT> const& q)
 {
-    const auto Ad = q.A() * r.dir;
-    const auto p = r.origin;
+    const auto Ad = q.A() * l.dir;
+    const auto p = l.pos;
 
     // Substituting x in Quadric equation x^TAx + 2b^Tx + c = 0 by ray equation x = t * dir + p yields
     // d^TAd t^2 + (2p^TAd + 2bd) t + p^TAp + 2bp + c = 0
-    const auto a = dot(r.dir, Ad);
-    const auto b = ScalarT(2) * (dot(p, Ad) + dot(q.b(), r.dir));
+    const auto a = dot(l.dir, Ad);
+    const auto b = ScalarT(2) * (dot(p, Ad) + dot(q.b(), l.dir));
     const auto c = dot(p, q.A() * vec3(p)) + ScalarT(2) * dot(q.b(), p) + q.c;
 
     // Solve the quadratic equation ax^2 + bx + c = 0
@@ -745,14 +776,7 @@ template <class ScalarT>
     const auto t2 = (-b + sqrtD) / (ScalarT(2) * a);
 
     const auto [tMin, tMax] = minmax(t1, t2);
-
-    if (tMin >= ScalarT(0))
-        return {tMin, tMax};
-
-    if (tMax >= ScalarT(0))
-        return tMax;
-
-    return {};
+    return {tMin, tMax};
 }
 
 
@@ -920,14 +944,6 @@ template <class ScalarT>
     return {};
 }
 
-template <class ScalarT>
-[[nodiscard]] constexpr ScalarT intersection_parameter(line<2, ScalarT> const& l0, line<2, ScalarT> const& l1)
-{
-    auto M = tg::mat<2, 2, ScalarT>::from_cols(l0.dir, -l1.dir);
-    auto t = inverse(M) * (l1.pos - l0.pos);
-    return t.x;
-}
-
 template <class ScalarT>
 [[nodiscard]] constexpr pair<ScalarT, ScalarT> intersection_parameters(line<2, ScalarT> const& l0, line<2, ScalarT> const& l1)
 {
@@ -960,12 +976,6 @@ template <int D, class ScalarT>
     return {res};
 }
 
-template <class ScalarT>
-[[nodiscard]] constexpr ScalarT intersection_parameter(line<3, ScalarT> const& l, plane<3, ScalarT> const& p)
-{
-    return (p.dis - dot(l.pos, p.normal)) / dot(l.dir, p.normal);
-}
-
 template <class ScalarT>
 [[nodiscard]] constexpr pos<3, ScalarT> intersection(plane<3, ScalarT> const& a, plane<3, ScalarT> const& b, plane<3, ScalarT> const& c)
 {
-- 
GitLab


From 7aa50e07ec78dc6a24c11222d17ed325dae94e55 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Wed, 18 Nov 2020 22:39:21 +0100
Subject: [PATCH 31/42] Renamed ray_hits to hits and ray_interval to
 hit_interval

---
 .../functions/objects/intersection.hh         | 84 +++++++++----------
 1 file changed, 42 insertions(+), 42 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index b25bcc4..4286689 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -32,7 +32,7 @@
 // intersects(a, b)              -> bool
 // intersection(a, b)            -> ???
 // intersection_safe(a, b)       -> optional<???>
-// intersection_parameter(a, b)  -> coords? (for a line or a ray: ray_hits<N, ScalarT> or optional<ray_interval> (when b is solid))
+// intersection_parameter(a, b)  -> coords? (for a line or a ray: hits<N, ScalarT> or optional<hit_interval> (when b is solid))
 // intersection_parameters(a, b) -> pair<coords, coords>?
 // intersection_exact(a, b)      -> variant
 // closest_intersection(a, b)            -> position/object (for a ray: optional<pos>)
@@ -80,14 +80,14 @@ namespace tg
 ///   operator[]
 ///   range-based-for
 template <int MaxHits, class HitT>
-struct ray_hits
+struct hits
 {
-    static constexpr bool is_ray_hits = true; // tag
+    static constexpr bool is_hits = true; // tag
     static constexpr int max_hits = MaxHits;
     using hit_t = HitT;
 
     template <class OtherT>
-    using as_ray_hits = ray_hits<MaxHits, OtherT>;
+    using as_hits = hits<MaxHits, OtherT>;
 
     [[nodiscard]] int size() const { return _size; }
     [[nodiscard]] bool any() const { return _size > 0; }
@@ -111,14 +111,14 @@ struct ray_hits
     [[nodiscard]] HitT const* begin() const { return _hit; }
     [[nodiscard]] HitT const* end() const { return _hit + _size; }
 
-    ray_hits() = default;
-    ray_hits(HitT* hits, int size) : _size(size)
+    hits() = default;
+    hits(HitT* hits, int size) : _size(size)
     {
         for (auto i = 0; i < size; ++i)
             _hit[i] = hits[i];
     }
     template <typename... HitTs>
-    ray_hits(HitTs... hits) : _size(sizeof...(HitTs)), _hit{hits...}
+    hits(HitTs... hits) : _size(sizeof...(HitTs)), _hit{hits...}
     {
     }
 
@@ -127,11 +127,11 @@ private:
     HitT _hit[MaxHits];
 };
 
-/// describes a continuous interval on a ray between start and end
+/// describes a continuous interval on a line or ray between start and end
 template <class ScalarT>
-struct ray_interval
+struct hit_interval
 {
-    static constexpr bool is_ray_interval = true; // tag
+    static constexpr bool is_hit_interval = true; // tag
 
     ScalarT start;
     ScalarT end;
@@ -146,7 +146,7 @@ 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>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> merge_hits(line<D, ScalarT> const& line, Objs const&... objs)
+[[nodiscard]] constexpr hits<2, ScalarT> merge_hits(line<D, ScalarT> const& line, Objs const&... objs)
 {
     ScalarT hits[2];
     hits[0] = tg::max<ScalarT>();
@@ -244,10 +244,10 @@ template <class A, class B>
     return {};
 }
 
-// if ray_hits intersection parameter is available, use that
+// if hits intersection parameter is available, use that
 template <class A, class B>
 [[nodiscard]] constexpr auto intersection(A const& a, B const& b) ->
-    typename decltype(intersection_parameter(a, b))::template as_ray_hits<typename A::pos_t>
+    typename decltype(intersection_parameter(a, b))::template as_hits<typename A::pos_t>
 {
     auto ts = intersection_parameter(a, b);
     typename A::pos_t hits[ts.max_hits];
@@ -256,10 +256,10 @@ template <class A, class B>
     return {hits, ts.size()};
 }
 
-// if an optional ray_interval intersection parameter is available, use that
+// if an optional hit_interval intersection parameter is available, use that
 template <class A, class B>
 [[nodiscard]] constexpr auto intersection(A const& a, B const& b)
-    -> std::enable_if<decltype(intersection_parameter(a, b).value())::is_ray_interval, // condition
+    -> std::enable_if<decltype(intersection_parameter(a, b).value())::is_hit_interval, // condition
     optional<segment<object_traits<A>::domain_dimension, typename object_traits<A>::scalar_t>>> // return type
 {
     auto ts = intersection_parameter(a, b);
@@ -268,7 +268,7 @@ template <class A, class B>
     return {};
 }
 
-// if ray_hits intersection parameter is available, use that
+// if hits intersection parameter is available, use that
 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>
 {
@@ -278,7 +278,7 @@ template <class A, class B>
     return {};
 }
 
-// if ray_interval intersection parameter is available, use that
+// if an optional hit_interval intersection parameter is available, use that
 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)>
 {
@@ -288,12 +288,12 @@ template <class A, class B>
     return {};
 }
 
-// if boundary_of a solid object returns ray_hits, use this to construct the ray_interval result of the solid intersection
+// if boundary_of a solid object returns hits, use this to construct the hit_interval result of the solid intersection
 template <int D, class ScalarT, class Obj>
 [[nodiscard]] constexpr auto intersection_parameter(line<D, ScalarT> const& l, Obj const& obj)
-    -> enable_if<!std::is_same_v<Obj, decltype(boundary_of(obj))>, optional<ray_interval<ScalarT>>>
+    -> enable_if<!std::is_same_v<Obj, decltype(boundary_of(obj))>, optional<hit_interval<ScalarT>>>
 {
-    const ray_hits<2, ScalarT> inter = intersection_parameter(l, boundary_of(obj));
+    const hits<2, ScalarT> inter = intersection_parameter(l, boundary_of(obj));
 
     if (inter.size() == 2)
         return {{inter[0], inter[1]}};
@@ -342,7 +342,7 @@ constexpr optional<pos<D, ScalarT>> intersection(Obj const& obj, pos<D, ScalarT>
 
 template <int D, class ScalarT, class Obj>
 [[nodiscard]] constexpr auto intersection_parameter(ray<D, ScalarT> const& ray, Obj const& obj)
-    -> decltype(intersection_parameter(inf_of(ray), obj).is_ray_hits, intersection_parameter(inf_of(ray), obj))
+    -> decltype(intersection_parameter(inf_of(ray), obj).is_hits, intersection_parameter(inf_of(ray), obj))
 {
     const auto inter = intersection_parameter(inf_of(ray), obj);
     constexpr auto maxHits = inter.max_hits;
@@ -361,7 +361,7 @@ template <int D, class ScalarT, class Obj>
 }
 template <int D, class ScalarT, class Obj>
 [[nodiscard]] constexpr auto intersection_parameter(ray<D, ScalarT> const& ray, Obj const& obj)
-    -> decltype(intersection_parameter(inf_of(ray), obj).value().is_ray_interval, intersection_parameter(inf_of(ray), obj))
+    -> decltype(intersection_parameter(inf_of(ray), obj).value().is_hit_interval, intersection_parameter(inf_of(ray), obj))
 {
     const auto inter = intersection_parameter(inf_of(ray), obj);
 
@@ -382,14 +382,14 @@ template <int D, class ScalarT, class Obj>
 
 // line - point
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<1, ScalarT> const& l, pos<1, ScalarT> const& p)
+[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<1, ScalarT> const& l, pos<1, ScalarT> const& p)
 {
     return coordinates(l, p);
 }
 
 // line - line
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<2, ScalarT> const& l0, line<2, ScalarT> const& l1)
+[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<2, ScalarT> const& l0, line<2, ScalarT> const& l1)
 {
     // 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);
@@ -399,7 +399,7 @@ template <class ScalarT>
 
 // line - ray
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<2, ScalarT> const& l, ray<2, ScalarT> const& r)
+[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<2, ScalarT> const& l, ray<2, ScalarT> const& r)
 {
     auto M = mat<2, 2, ScalarT>::from_cols(l.dir, -r.dir);
     auto t = inverse(M) * (r.origin - l.pos);
@@ -410,7 +410,7 @@ template <class ScalarT>
 
 // line - segment
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<2, ScalarT> const& l, segment<2, ScalarT> const& s)
+[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<2, ScalarT> const& l, segment<2, ScalarT> const& s)
 {
     auto M = mat<2, 2, ScalarT>::from_cols(l.dir, s.pos0 - s.pos1);
     auto t = inverse(M) * (s.pos0 - l.pos);
@@ -421,7 +421,7 @@ template <class ScalarT>
 
 // line - plane
 template <int D, class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<D, ScalarT> const& l, plane<D, ScalarT> const& p)
+[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<D, ScalarT> const& l, plane<D, ScalarT> const& p)
 {
     const auto dotND = dot(p.normal, l.dir);
     if (dotND == ScalarT(0)) // if plane normal and line direction are orthogonal, there is no intersection
@@ -433,7 +433,7 @@ template <int D, class ScalarT>
 
 // line - halfspace
 template <int D, class ScalarT>
-[[nodiscard]] constexpr optional<ray_interval<ScalarT>> intersection_parameter(line<D, ScalarT> const& l, halfspace<D, ScalarT> const& h)
+[[nodiscard]] constexpr optional<hit_interval<ScalarT>> intersection_parameter(line<D, ScalarT> const& l, halfspace<D, ScalarT> const& h)
 {
     const auto dotND = dot(h.normal, l.dir);
     const auto dist = signed_distance(l.pos, h);
@@ -468,7 +468,7 @@ template <int D, class ScalarT>
 
 // line - aabb
 template <int D, class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, aabb_boundary<D, ScalarT> const& b)
+[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, aabb_boundary<D, ScalarT> const& b)
 {
     auto tFirst = tg::min<ScalarT>();
     auto tSecond = tg::max<ScalarT>();
@@ -494,7 +494,7 @@ template <int D, class ScalarT>
 
 // line - box
 template <int D, class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, box_boundary<D, ScalarT> const& b)
+[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, box_boundary<D, ScalarT> const& b)
 {
     const auto bMin = b[comp<D, ScalarT>(-1)] - l.pos;
     const auto bMax = b[comp<D, ScalarT>(1)] - l.pos;
@@ -521,7 +521,7 @@ template <int D, class ScalarT>
     return {tFirst, tSecond};
 }
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<3, ScalarT> const& l, box<2, ScalarT, 3> const& b)
+[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<3, ScalarT> const& l, box<2, ScalarT, 3> const& b)
 {
     const auto t = intersection_parameter(l, plane<3, ScalarT>(normal_of(b), b.center));
     if (!t.any()) // line parallel to box plane
@@ -536,7 +536,7 @@ template <class ScalarT>
 
 // line - disk
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<3, ScalarT> const& l, sphere<2, ScalarT, 3> const& d)
+[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<3, ScalarT> const& l, sphere<2, ScalarT, 3> const& d)
 {
     const auto t = intersection_parameter(l, plane<3, ScalarT>(d.normal, d.center));
     if (!t.any()) // line parallel to disk plane
@@ -551,7 +551,7 @@ template <class ScalarT>
 
 // line - sphere
 template <int D, class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, sphere_boundary<D, ScalarT> const& s)
+[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, sphere_boundary<D, ScalarT> const& s)
 {
     auto t = dot(s.center - l.pos, l.dir);
 
@@ -566,12 +566,12 @@ template <int D, class ScalarT>
 
 // line - hemisphere
 template <int D, class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, hemisphere_boundary<D, ScalarT> const& h)
+[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, hemisphere_boundary<D, ScalarT> const& h)
 {
     return detail::merge_hits(l, caps_of(h), boundary_no_caps_of(h));
 }
 template <int D, class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, hemisphere_boundary_no_caps<D, ScalarT> const& h)
+[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, hemisphere_boundary_no_caps<D, ScalarT> const& h)
 {
     ScalarT hits[2];
     auto numHits = 0;
@@ -595,7 +595,7 @@ template <class Obj, int D, class ScalarT, class TraitsT>
 
 // line - capsule
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, capsule_boundary<3, ScalarT> const& c)
+[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, capsule_boundary<3, ScalarT> const& c)
 {
     using caps_t = hemisphere_boundary_no_caps<3, ScalarT>;
     const auto n = direction(c);
@@ -613,7 +613,7 @@ template <class Obj, class ScalarT, class TraitsT>
 
 // line - cylinder
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, cylinder_boundary<3, ScalarT> const& c)
+[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, cylinder_boundary<3, ScalarT> const& c)
 {
     const auto caps = caps_of(c);
     return detail::merge_hits(l, caps[0], caps[1], boundary_no_caps_of(c));
@@ -631,7 +631,7 @@ template <class Obj, class ScalarT, class TraitsT>
     }
 }
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, cylinder_boundary_no_caps<3, ScalarT> const& c)
+[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, cylinder_boundary_no_caps<3, ScalarT> const& c)
 {
     auto cdir = direction(c);
     auto cosA = dot(cdir, l.dir);
@@ -681,7 +681,7 @@ template <class ScalarT>
 
 // line - inf_cylinder
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<2, ScalarT> const& l, inf_cylinder_boundary<2, ScalarT> const& c)
+[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<2, ScalarT> const& l, inf_cylinder_boundary<2, ScalarT> const& c)
 {
     const auto n = perpendicular(c.axis.dir);
     const auto d = dot(l.dir, n);
@@ -695,7 +695,7 @@ template <class ScalarT>
 
 // line - triangle2
 template <class ScalarT>
-[[nodiscard]] constexpr optional<ray_interval<ScalarT>> intersection_parameter(line<2, ScalarT> const& l, triangle<2, ScalarT> const& t)
+[[nodiscard]] constexpr optional<hit_interval<ScalarT>> intersection_parameter(line<2, ScalarT> const& l, triangle<2, ScalarT> const& t)
 {
     ScalarT closestIntersection = tg::max<ScalarT>();
     ScalarT furtherIntersection = tg::min<ScalarT>();
@@ -720,7 +720,7 @@ template <class ScalarT>
 
 // line - triangle3
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<1, ScalarT> intersection_parameter(line<3, ScalarT> const& l,
+[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<3, ScalarT> const& l,
                                                                  triangle<3, ScalarT> const& t,
                                                                  dont_deduce<ScalarT> eps = 100 * tg::epsilon<ScalarT>)
 {
@@ -755,7 +755,7 @@ template <class ScalarT>
 
 // line - quadric_boundary (as an isosurface, not error quadric)
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, quadric<3, ScalarT> const& q)
+[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, quadric<3, ScalarT> const& q)
 {
     const auto Ad = q.A() * l.dir;
     const auto p = l.pos;
-- 
GitLab


From b2642eb7a1630106a71bf50e46ba6599830e226b Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Wed, 18 Nov 2020 23:04:48 +0100
Subject: [PATCH 32/42] Fixed gcc error

---
 src/typed-geometry/types/objects/line.hh    | 4 ++--
 src/typed-geometry/types/objects/ray.hh     | 4 ++--
 src/typed-geometry/types/objects/segment.hh | 1 +
 3 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/typed-geometry/types/objects/line.hh b/src/typed-geometry/types/objects/line.hh
index 02d2891..80071c7 100644
--- a/src/typed-geometry/types/objects/line.hh
+++ b/src/typed-geometry/types/objects/line.hh
@@ -46,8 +46,8 @@ template <int D, class ScalarT>
 struct line
 {
     using scalar_t = ScalarT;
-    using dir_t = dir<D, ScalarT>;
-    using pos_t = pos<D, ScalarT>;
+    using dir_t = tg::dir<D, ScalarT>; // Note: gcc requires the tg:: prefix to distinguish from member variable
+    using pos_t = tg::pos<D, ScalarT>;
 
     pos_t pos;
     dir_t dir;
diff --git a/src/typed-geometry/types/objects/ray.hh b/src/typed-geometry/types/objects/ray.hh
index a376315..f0b85df 100644
--- a/src/typed-geometry/types/objects/ray.hh
+++ b/src/typed-geometry/types/objects/ray.hh
@@ -44,8 +44,8 @@ template <int D, class ScalarT>
 struct ray
 {
     using scalar_t = ScalarT;
-    using dir_t = dir<D, ScalarT>;
-    using pos_t = pos<D, ScalarT>;
+    using dir_t = tg::dir<D, ScalarT>; // Note: gcc requires the tg:: prefix to distinguish from member variable
+    using pos_t = tg::pos<D, ScalarT>;
 
     pos_t origin;
     dir_t dir;
diff --git a/src/typed-geometry/types/objects/segment.hh b/src/typed-geometry/types/objects/segment.hh
index 8f6f91a..d80e584 100644
--- a/src/typed-geometry/types/objects/segment.hh
+++ b/src/typed-geometry/types/objects/segment.hh
@@ -44,6 +44,7 @@ using usegment4 = segment<4, u32>;
 template <int D, class ScalarT>
 struct segment
 {
+    using scalar_t = ScalarT;
     using pos_t = pos<D, ScalarT>;
 
     pos_t pos0;
-- 
GitLab


From 3fc80be390fde45195ec165e4cb4635d069a6943 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Thu, 19 Nov 2020 15:53:15 +0100
Subject: [PATCH 33/42] Added inf_cylinder intersection

---
 .../functions/objects/direction.hh            |  5 ++
 .../functions/objects/intersection.hh         | 78 +++++++++----------
 2 files changed, 44 insertions(+), 39 deletions(-)

diff --git a/src/typed-geometry/functions/objects/direction.hh b/src/typed-geometry/functions/objects/direction.hh
index 858676a..79c3d98 100644
--- a/src/typed-geometry/functions/objects/direction.hh
+++ b/src/typed-geometry/functions/objects/direction.hh
@@ -44,6 +44,11 @@ template <int D, class ScalarT, class TraitsT>
     return direction(c.axis);
 }
 template <int D, class ScalarT, class TraitsT>
+[[nodiscard]] constexpr dir<D, ScalarT> direction(inf_cylinder<D, ScalarT, TraitsT> const& c)
+{
+    return direction(c.axis);
+}
+template <int D, class ScalarT, class TraitsT>
 [[nodiscard]] constexpr dir<D, ScalarT> direction(capsule<D, ScalarT, TraitsT> const& c)
 {
     return direction(c.axis);
diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 4286689..f2c8466 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -176,7 +176,7 @@ template <class Obj, class... Objs>
 }
 
 template <class A, class B>
-using try_closest_intersection_parameter = decltype(closest_intersection_parameter(std::declval<A const&>(), std::declval<B const&>()).has_value());
+using try_closest_intersection_parameter = decltype(closest_intersection_parameter(std::declval<A const&>(), std::declval<B const&>()));
 }
 
 
@@ -633,53 +633,53 @@ template <class Obj, class ScalarT, class TraitsT>
 template <class ScalarT>
 [[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, cylinder_boundary_no_caps<3, ScalarT> const& c)
 {
-    auto cdir = direction(c);
-    auto cosA = dot(cdir, l.dir);
-    auto sinA_sqr = 1 - cosA * cosA;
+    auto const infInter = intersection_parameter(l, inf_of(c));
+    if (!infInter.any())
+        return infInter;
 
-    if (sinA_sqr <= 0)
-        return {}; // line and cylinder are parallel
+    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);
 
-    // compute closest points of line and cylinder axis
-    auto origDiff = l.pos - c.axis.pos0;
-    auto fLine = dot(l.dir, origDiff);
-    auto fAxis = dot(cdir, origDiff);
-    auto tLine = (cosA * fAxis - fLine) / sinA_sqr;
-    auto tAxis = (fAxis - cosA * fLine) / sinA_sqr;
-
-    auto closest_on_line = l.pos + tLine * l.dir;
-    auto closest_on_axis = c.axis.pos0 + tAxis * cdir;
-    auto line_axis_dist_sqr = distance_sqr(closest_on_line, closest_on_axis);
-    auto cyl_radius_sqr = c.radius * c.radius;
-
-    if (line_axis_dist_sqr > cyl_radius_sqr)
-        return {}; // line misses the cylinder
+    ScalarT hits[2];
+    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];
 
-    // radius in 2D slice
-    auto r_2D = sqrt(cyl_radius_sqr - line_axis_dist_sqr);
+    return {hits, numHits};
+}
 
-    // infinite tube intersection
-    auto s = r_2D / sqrt(sinA_sqr);
-    auto cyl_intersection_0 = closest_on_line - s * l.dir;
-    auto cyl_intersection_1 = closest_on_line + s * l.dir;
+// 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;
 
-    // project onto line segment
-    auto line_length = length(c.axis);
-    auto lambda_0 = dot(cyl_intersection_0 - c.axis.pos0, cdir);
-    auto lambda_1 = dot(cyl_intersection_1 - c.axis.pos0, cdir);
+    if (sinA_sqr <= 0)
+        return {}; // line and cylinder are parallel
 
-    ScalarT hits[2];
-    int hit_cnt = 0;
+    // 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;
 
-    if (0 <= lambda_0 && lambda_0 < line_length)
-        hits[hit_cnt++] = tLine - s;
-    if (0 <= lambda_1 && lambda_1 < line_length)
-        hits[hit_cnt++] = tLine + s;
+    auto const line_axis_dist_sqr = distance_sqr(l[tLine], c.axis[tAxis]);
+    auto const cyl_radius_sqr = c.radius * c.radius;
 
-    return {hits, hit_cnt};
-}
+    if (cyl_radius_sqr < line_axis_dist_sqr)
+        return {}; // line misses the cylinder
 
-// line - inf_cylinder
+    // 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};
+}
 template <class ScalarT>
 [[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<2, ScalarT> const& l, inf_cylinder_boundary<2, ScalarT> const& c)
 {
-- 
GitLab


From 3742a0814ba9b6f80f762ba480ac96b2c6374386 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Thu, 19 Nov 2020 16:19:11 +0100
Subject: [PATCH 34/42] Added inf_cone2 and inf_cylinder2 common types

---
 src/typed-geometry/types/objects/inf_cone.hh     | 6 ++++++
 src/typed-geometry/types/objects/inf_cylinder.hh | 6 ++++++
 2 files changed, 12 insertions(+)

diff --git a/src/typed-geometry/types/objects/inf_cone.hh b/src/typed-geometry/types/objects/inf_cone.hh
index 6e17163..1af69f0 100644
--- a/src/typed-geometry/types/objects/inf_cone.hh
+++ b/src/typed-geometry/types/objects/inf_cone.hh
@@ -16,6 +16,12 @@ template <int D, class ScalarT, class TraitsT = default_object_tag>
 struct inf_cone;
 
 // Common cone types
+using inf_cone2 = inf_cone<2, f32>;
+using finf_cone2 = inf_cone<2, f32>;
+using dinf_cone2 = inf_cone<2, f64>;
+using iinf_cone2 = inf_cone<2, i32>;
+using uinf_cone2 = inf_cone<2, u32>;
+
 using inf_cone3 = inf_cone<3, f32>;
 using finf_cone3 = inf_cone<3, f32>;
 using dinf_cone3 = inf_cone<3, f64>;
diff --git a/src/typed-geometry/types/objects/inf_cylinder.hh b/src/typed-geometry/types/objects/inf_cylinder.hh
index f4c1a16..2590022 100644
--- a/src/typed-geometry/types/objects/inf_cylinder.hh
+++ b/src/typed-geometry/types/objects/inf_cylinder.hh
@@ -16,6 +16,12 @@ template <int D, class ScalarT, class TraitsT = default_object_tag>
 struct inf_cylinder;
 
 // Common infinite cylinder types
+using inf_cylinder2 = inf_cylinder<2, f32>;
+using finf_cylinder2 = inf_cylinder<2, f32>;
+using dinf_cylinder2 = inf_cylinder<2, f64>;
+using iinf_cylinder2 = inf_cylinder<2, i32>;
+using uinf_cylinder2 = inf_cylinder<2, u32>;
+
 using inf_cylinder3 = inf_cylinder<3, f32>;
 using finf_cylinder3 = inf_cylinder<3, f32>;
 using dinf_cylinder3 = inf_cylinder<3, f64>;
-- 
GitLab


From a59574fb8a7d35211128c78629d265b02bc90c96 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Thu, 19 Nov 2020 23:03:42 +0100
Subject: [PATCH 35/42] Added inf_cone intersections

---
 .../functions/objects/any_point.hh            |  7 +-
 .../functions/objects/intersection.hh         | 64 +++++++++++++++----
 2 files changed, 58 insertions(+), 13 deletions(-)

diff --git a/src/typed-geometry/functions/objects/any_point.hh b/src/typed-geometry/functions/objects/any_point.hh
index bdf782d..43a27b3 100644
--- a/src/typed-geometry/functions/objects/any_point.hh
+++ b/src/typed-geometry/functions/objects/any_point.hh
@@ -98,7 +98,12 @@ template <int D, class ScalarT, class TraitsT>
     if constexpr (std::is_same_v<TraitsT, default_object_tag>)
         return c.axis.pos;
     else
-        return c.axis.pos + c.radius * any_normal(c.axis.dir);
+    {
+        if constexpr (D == 3)
+            return c.axis.pos + c.radius * any_normal(c.axis.dir);
+        else
+            return c.axis.pos + c.radius * perpendicular(c.axis.dir);
+    }
 }
 
 template <int D, class ScalarT>
diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index f2c8466..bf94a66 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -10,6 +10,8 @@
 #include <typed-geometry/types/objects/cylinder.hh>
 #include <typed-geometry/types/objects/halfspace.hh>
 #include <typed-geometry/types/objects/hemisphere.hh>
+#include <typed-geometry/types/objects/inf_cone.hh>
+#include <typed-geometry/types/objects/inf_cylinder.hh>
 #include <typed-geometry/types/objects/line.hh>
 #include <typed-geometry/types/objects/plane.hh>
 #include <typed-geometry/types/objects/ray.hh>
@@ -175,6 +177,22 @@ template <class Obj, class... Objs>
     return (intersects(obj, objs) || ...);
 }
 
+// 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};
+}
+
 template <class A, class B>
 using try_closest_intersection_parameter = decltype(closest_intersection_parameter(std::declval<A const&>(), std::declval<B const&>()));
 }
@@ -470,6 +488,7 @@ template <int D, class ScalarT>
 template <int D, class ScalarT>
 [[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, aabb_boundary<D, ScalarT> const& b)
 {
+    // based on ideas from https://gamedev.stackexchange.com/q/18436
     auto tFirst = tg::min<ScalarT>();
     auto tSecond = tg::max<ScalarT>();
     for (auto i = 0; i < D; ++i)
@@ -693,6 +712,38 @@ template <class ScalarT>
     return {tMin, tMax};
 }
 
+// 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);
+
+    // exclude intersections with mirrored cone
+    ScalarT hits[2];
+    auto numHits = 0;
+    if (dot(l[inter[0]] - ic.apex, ic.opening_dir) >= ScalarT(0))
+        hits[numHits++] = inter[0];
+    if (dot(l[inter[1]] - ic.apex, ic.opening_dir) >= ScalarT(0))
+        hits[numHits++] = inter[1];
+
+    return {hits, numHits};
+}
+
 // line - triangle2
 template <class ScalarT>
 [[nodiscard]] constexpr optional<hit_interval<ScalarT>> intersection_parameter(line<2, ScalarT> const& l, triangle<2, ScalarT> const& t)
@@ -765,18 +816,7 @@ template <class ScalarT>
     const auto a = dot(l.dir, Ad);
     const auto b = ScalarT(2) * (dot(p, Ad) + dot(q.b(), l.dir));
     const auto c = dot(p, q.A() * vec3(p)) + ScalarT(2) * dot(q.b(), p) + q.c;
-
-    // Solve the quadratic equation ax^2 + bx + c = 0
-    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};
+    return detail::solve_quadratic(a, b, c);
 }
 
 
-- 
GitLab


From 05a25760acd58e2a87d3627a2f6d3d9f40ecc95e Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Fri, 20 Nov 2020 02:17:01 +0100
Subject: [PATCH 36/42] Added ellipse2in3 intersection

---
 .../functions/objects/intersection.hh         | 28 +++++++++++++++++--
 .../functions/objects/normal.hh               |  7 +++++
 2 files changed, 33 insertions(+), 2 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index bf94a66..cf18b93 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -8,6 +8,7 @@
 #include <typed-geometry/types/objects/box.hh>
 #include <typed-geometry/types/objects/capsule.hh>
 #include <typed-geometry/types/objects/cylinder.hh>
+#include <typed-geometry/types/objects/ellipse.hh>
 #include <typed-geometry/types/objects/halfspace.hh>
 #include <typed-geometry/types/objects/hemisphere.hh>
 #include <typed-geometry/types/objects/inf_cone.hh>
@@ -732,13 +733,16 @@ template <class ScalarT>
     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
     ScalarT hits[2];
     auto numHits = 0;
-    if (dot(l[inter[0]] - ic.apex, ic.opening_dir) >= ScalarT(0))
+    auto const coneDir = ic.opening_angle > 180_deg ? -ic.opening_dir : ic.opening_dir;
+    if (dot(l[inter[0]] - ic.apex, coneDir) >= ScalarT(0))
         hits[numHits++] = inter[0];
-    if (dot(l[inter[1]] - ic.apex, ic.opening_dir) >= ScalarT(0))
+    if (dot(l[inter[1]] - ic.apex, coneDir) >= ScalarT(0))
         hits[numHits++] = inter[1];
 
     return {hits, numHits};
@@ -804,6 +808,26 @@ template <class ScalarT>
     return (ScalarT(1) / det) * dot(e2, qvec);
 }
 
+// line - ellipse
+template <class ScalarT>
+[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<3, ScalarT> const& l, ellipse<2, ScalarT, 3> const& e)
+{
+    const auto t = intersection_parameter(l, plane<3, ScalarT>(normal_of(e), e.center));
+    if (!t.any()) // line parallel to ellipse plane
+        return {};
+
+    // simplified contains(e, p) without plane check and eps == 0
+    auto pc = l[t.first()] - e.center;
+    auto x = dot(pc, e.semi_axes[0]);
+    auto y = dot(pc, e.semi_axes[1]);
+    auto a = length_sqr(e.semi_axes[0]);
+    auto b = length_sqr(e.semi_axes[1]);
+
+    if (pow2(x / a) + pow2(y / b) <= ScalarT(1))
+        return t;
+    return {};
+}
+
 // line - quadric_boundary (as an isosurface, not error quadric)
 template <class ScalarT>
 [[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, quadric<3, ScalarT> const& q)
diff --git a/src/typed-geometry/functions/objects/normal.hh b/src/typed-geometry/functions/objects/normal.hh
index dd09e8b..9096e06 100644
--- a/src/typed-geometry/functions/objects/normal.hh
+++ b/src/typed-geometry/functions/objects/normal.hh
@@ -1,6 +1,7 @@
 #pragma once
 
 #include <typed-geometry/types/objects/box.hh>
+#include <typed-geometry/types/objects/ellipse.hh>
 #include <typed-geometry/types/objects/halfspace.hh>
 #include <typed-geometry/types/objects/line.hh>
 #include <typed-geometry/types/objects/plane.hh>
@@ -78,6 +79,12 @@ template <class ScalarT, class TraitsT>
     return normalize(cross(b.half_extents[0], b.half_extents[1]));
 }
 
+template <class ScalarT, class TraitsT>
+[[nodiscard]] constexpr dir<3, fractional_result<ScalarT>> normal_of(ellipse<2, ScalarT, 3, TraitsT> const& e)
+{
+    return normalize(cross(e.semi_axes[0], e.semi_axes[1]));
+}
+
 template <class ScalarT>
 [[nodiscard]] constexpr dir<3, ScalarT> any_normal(vec<3, ScalarT> const& v)
 {
-- 
GitLab


From ef0d532b8b9c689bfa366d7da42e5ed1593be15f Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Fri, 20 Nov 2020 13:20:15 +0100
Subject: [PATCH 37/42] Added ellipse intersections

---
 .../functions/objects/intersection.hh         | 25 ++++++++++++++++---
 1 file changed, 22 insertions(+), 3 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index cf18b93..b48459a 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -739,10 +739,13 @@ template <class ScalarT>
     // exclude intersections with mirrored cone
     ScalarT hits[2];
     auto numHits = 0;
-    auto const coneDir = ic.opening_angle > 180_deg ? -ic.opening_dir : ic.opening_dir;
-    if (dot(l[inter[0]] - ic.apex, coneDir) >= ScalarT(0))
+    TG_ASSERT(ic.opening_angle <= 180_deg, "Only convex objects are supported, but an inf_cone with openinge angle > 180 degree is not convex.");
+    // 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))
         hits[numHits++] = inter[0];
-    if (dot(l[inter[1]] - ic.apex, coneDir) >= ScalarT(0))
+    if (dot(l[inter[1]] - ic.apex, ic.opening_dir) >= ScalarT(0))
         hits[numHits++] = inter[1];
 
     return {hits, numHits};
@@ -827,6 +830,22 @@ template <class ScalarT>
         return t;
     return {};
 }
+template <int D, class ScalarT>
+[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, ellipse_boundary<D, ScalarT> const& e)
+{
+    // transform line to ellipse space (ellipse gets unit sphere at origin)
+    auto const pc = l.pos - e.center;
+    pos<D, ScalarT> p;
+    vec<D, ScalarT> d; // in ellipse space, this is no longer a unit vector
+    for (auto i = 0; i < D; ++i)
+    {
+        auto const axis2 = dot(e.semi_axes[i], e.semi_axes[i]);
+        p[i] = dot(pc, e.semi_axes[i]) / axis2;
+        d[i] = dot(l.dir, e.semi_axes[i]) / axis2;
+    }
+    // to find intersection with unit sphere, <t*d + p, t*d + p> == 1 has to be solved
+    return detail::solve_quadratic(dot(d, d), ScalarT(2) * dot(d, p), dot(p, p) - ScalarT(1));
+}
 
 // line - quadric_boundary (as an isosurface, not error quadric)
 template <class ScalarT>
-- 
GitLab


From 674fb112732a1a08d5abd237abc0f2a45580ad7f Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Mon, 7 Dec 2020 22:39:34 +0100
Subject: [PATCH 38/42] Added intersections for all pyramids

---
 .../functions/objects/intersection.hh         | 56 ++++++++++++++++++-
 1 file changed, 55 insertions(+), 1 deletion(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index b48459a..7bd3ea8 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -15,6 +15,7 @@
 #include <typed-geometry/types/objects/inf_cylinder.hh>
 #include <typed-geometry/types/objects/line.hh>
 #include <typed-geometry/types/objects/plane.hh>
+#include <typed-geometry/types/objects/pyramid.hh>
 #include <typed-geometry/types/objects/ray.hh>
 #include <typed-geometry/types/objects/segment.hh>
 #include <typed-geometry/types/objects/sphere.hh>
@@ -28,6 +29,7 @@
 #include "contains.hh"
 #include "coordinates.hh"
 #include "direction.hh"
+#include "faces.hh"
 #include "normal.hh"
 
 // family of intersection functions:
@@ -170,6 +172,11 @@ template <int D, class ScalarT, class... Objs>
     TG_ASSERT(numHits <= 2);
     return {hits, numHits};
 }
+template <int D, class ScalarT, class ObjT, int N, std::size_t... I>
+[[nodiscard]] constexpr hits<2, ScalarT> merge_hits(line<D, ScalarT> const& line, array<ObjT, N> objs, std::index_sequence<I...>)
+{
+    return merge_hits(line, objs[I]...);
+}
 
 // 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>
@@ -739,7 +746,7 @@ template <class ScalarT>
     // exclude intersections with mirrored cone
     ScalarT hits[2];
     auto numHits = 0;
-    TG_ASSERT(ic.opening_angle <= 180_deg, "Only convex objects are supported, but an inf_cone with openinge angle > 180 degree is not convex.");
+    TG_ASSERT(ic.opening_angle <= 180_deg && "Only convex objects are supported, but an inf_cone with openinge angle > 180 degree is not convex.");
     // 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)) ...
@@ -751,6 +758,53 @@ template <class ScalarT>
     return {hits, numHits};
 }
 
+// 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
+    ScalarT hits[2];
+    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);
+    return detail::merge_hits(l, faces, std::make_index_sequence<faces.size()>{});
+}
+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));
+}
+
 // line - triangle2
 template <class ScalarT>
 [[nodiscard]] constexpr optional<hit_interval<ScalarT>> intersection_parameter(line<2, ScalarT> const& l, triangle<2, ScalarT> const& t)
-- 
GitLab


From b1cf5f84254d2ea8174fb6108ab5318fbf17b5aa Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Tue, 8 Dec 2020 01:29:01 +0100
Subject: [PATCH 39/42] Included <utility> header

---
 src/typed-geometry/functions/objects/intersection.hh | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 7bd3ea8..5e51cb2 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -32,6 +32,8 @@
 #include "faces.hh"
 #include "normal.hh"
 
+#include <utility>
+
 // family of intersection functions:
 
 // intersects(a, b)              -> bool
-- 
GitLab


From 7cd1c77016ccc510d685b57e3b3222405a73c1a6 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Tue, 8 Dec 2020 01:55:29 +0100
Subject: [PATCH 40/42] Removed deprecated usage of dot(pos, pos)

---
 src/typed-geometry/functions/objects/intersection.hh | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 5e51cb2..e870c17 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -891,7 +891,7 @@ template <int D, class ScalarT>
 {
     // transform line to ellipse space (ellipse gets unit sphere at origin)
     auto const pc = l.pos - e.center;
-    pos<D, ScalarT> p;
+    vec<D, ScalarT> p;
     vec<D, ScalarT> d; // in ellipse space, this is no longer a unit vector
     for (auto i = 0; i < D; ++i)
     {
-- 
GitLab


From a08acaeb56596de998b6b55d557e90f7cc71294d Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Tue, 8 Dec 2020 02:05:55 +0100
Subject: [PATCH 41/42] Fixed CI

---
 src/typed-geometry/functions/objects/intersection.hh | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index e870c17..2424529 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -174,8 +174,8 @@ template <int D, class ScalarT, class... Objs>
     TG_ASSERT(numHits <= 2);
     return {hits, numHits};
 }
-template <int D, class ScalarT, class ObjT, int N, std::size_t... I>
-[[nodiscard]] constexpr hits<2, ScalarT> merge_hits(line<D, ScalarT> const& line, array<ObjT, N> objs, std::index_sequence<I...>)
+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...>)
 {
     return merge_hits(line, objs[I]...);
 }
@@ -748,7 +748,7 @@ template <class ScalarT>
     // exclude intersections with mirrored cone
     ScalarT hits[2];
     auto numHits = 0;
-    TG_ASSERT(ic.opening_angle <= 180_deg && "Only convex objects are supported, but an inf_cone with openinge angle > 180 degree is not convex.");
+    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.");
     // 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)) ...
@@ -799,7 +799,7 @@ template <class BaseT, typename = std::enable_if_t<!std::is_same_v<BaseT, sphere
 [[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);
-    return detail::merge_hits(l, faces, std::make_index_sequence<faces.size()>{});
+    return detail::merge_hits_array(l, faces, std::make_index_sequence<faces.size()>{});
 }
 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)
-- 
GitLab


From 1a1a37b244af85f4fc28bdbb964b1065f3189d86 Mon Sep 17 00:00:00 2001
From: Aaron Grabowy <aaron.grabowy@rwth-aachen.de>
Date: Tue, 8 Dec 2020 11:34:22 +0100
Subject: [PATCH 42/42] Trying to fix gcc warning

---
 .../functions/objects/intersection.hh              | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/typed-geometry/functions/objects/intersection.hh b/src/typed-geometry/functions/objects/intersection.hh
index 2424529..709497e 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -131,7 +131,7 @@ struct hits
 
 private:
     int _size = 0;
-    HitT _hit[MaxHits];
+    HitT _hit[MaxHits] = {};
 };
 
 /// describes a continuous interval on a line or ray between start and end
@@ -155,7 +155,7 @@ namespace detail
 template <int D, class ScalarT, class... Objs>
 [[nodiscard]] constexpr hits<2, ScalarT> merge_hits(line<D, ScalarT> const& line, Objs const&... objs)
 {
-    ScalarT hits[2];
+    ScalarT hits[2] = {};
     hits[0] = tg::max<ScalarT>();
     hits[1] = tg::min<ScalarT>();
     auto numHits = 0;
@@ -278,7 +278,7 @@ template <class A, class B>
     typename decltype(intersection_parameter(a, b))::template as_hits<typename A::pos_t>
 {
     auto ts = intersection_parameter(a, b);
-    typename A::pos_t hits[ts.max_hits];
+    typename A::pos_t hits[ts.max_hits] = {};
     for (auto i = 0; i < ts.size(); ++i)
         hits[i] = a[ts[i]];
     return {hits, ts.size()};
@@ -602,7 +602,7 @@ template <int D, class ScalarT>
 template <int D, class ScalarT>
 [[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<D, ScalarT> const& l, hemisphere_boundary_no_caps<D, ScalarT> const& h)
 {
-    ScalarT hits[2];
+    ScalarT hits[2] = {};
     auto numHits = 0;
     const auto sphereHits = intersection_parameter(l, sphere_boundary<D, ScalarT>(h.center, h.radius));
     const auto halfSpace = halfspace<D, ScalarT>(-h.normal, h.center); // the intersection of this halfspace and the sphere is exactly the hemisphere
@@ -670,7 +670,7 @@ template <class ScalarT>
     auto const lambda0 = dot(l[infInter[0]] - c.axis.pos0, d);
     auto const lambda1 = dot(l[infInter[1]] - c.axis.pos0, d);
 
-    ScalarT hits[2];
+    ScalarT hits[2] = {};
     auto numHits = 0;
     auto const dDotD = dot(d, d);
     if (ScalarT(0) <= lambda0 && lambda0 <= dDotD)
@@ -746,7 +746,7 @@ template <class ScalarT>
         return inter;
 
     // exclude intersections with mirrored cone
-    ScalarT hits[2];
+    ScalarT hits[2] = {};
     auto numHits = 0;
     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.");
     // if it is not used for solid cones, this works:
@@ -782,7 +782,7 @@ template <class ScalarT>
         return inter;
 
     // exclude intersections with mirrored cone
-    ScalarT hits[2];
+    ScalarT hits[2] = {};
     auto numHits = 0;
     auto const h0 = dot(l[inter[0]] - apex, openingDir);
     auto const h1 = dot(l[inter[1]] - apex, openingDir);
-- 
GitLab