diff --git a/src/typed-geometry/detail/optional.hh b/src/typed-geometry/detail/optional.hh
index f7f3fd7e7a411aa690be267399eb2d8a19978df3..dba24b1e5b21e302ec51a591b40e0e30a573f541 100644
--- a/src/typed-geometry/detail/optional.hh
+++ b/src/typed-geometry/detail/optional.hh
@@ -45,7 +45,7 @@ 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)
 {
     return o.has_value() ? (os << o.value()) : (os << "[empty]");
 }
diff --git a/src/typed-geometry/detail/utility.hh b/src/typed-geometry/detail/utility.hh
index 4e5012daf15352c187710c00ddb4587a0000324c..d084e15053b990e8cb07e0e44a2e321e49ecb500 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/basic/scalar_math.hh b/src/typed-geometry/functions/basic/scalar_math.hh
index 59f4d243838313acdbd32c5066cab425f2a3d64a..ef841c5fe4804c401402ffe4b0f08613e47a78ea 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/any_point.hh b/src/typed-geometry/functions/objects/any_point.hh
index bdf782d4f07f62b4227b5d96e2a6041de74e8e05..43a27b3d8709015634a169e58d29a908445993cb 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/contains.hh b/src/typed-geometry/functions/objects/contains.hh
index 123d9efdd666cf0d8abb4baf0080196215d97934..25cecc8ec7e3ae87010fd9e6eae4b13548349abb 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))
 {
@@ -260,7 +274,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;
 }
diff --git a/src/typed-geometry/functions/objects/direction.hh b/src/typed-geometry/functions/objects/direction.hh
index 858676a213086184b7d1e3b861acd61ddaa30264..79c3d98f6809aabccf722acb5b319af08d0b27f2 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 98425f5c3874782e686361ef224e68c5166f4fbd..709497ef44950418dd301f1eac1280222a446c3e 100644
--- a/src/typed-geometry/functions/objects/intersection.hh
+++ b/src/typed-geometry/functions/objects/intersection.hh
@@ -4,9 +4,18 @@
 #include <typed-geometry/feature/assert.hh>
 #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/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>
+#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>
@@ -17,26 +26,29 @@
 #include <typed-geometry/functions/vector/length.hh>
 
 #include "aabb.hh"
-#include "closest_points.hh"
 #include "contains.hh"
+#include "coordinates.hh"
 #include "direction.hh"
+#include "faces.hh"
 #include "normal.hh"
 
+#include <utility>
+
 // family of intersection functions:
 
 // intersects(a, b)              -> bool
 // intersection(a, b)            -> ???
 // intersection_safe(a, b)       -> optional<???>
-// intersection_parameter(a, b)  -> coords?
+// 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)            -> ???
-// 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 +58,24 @@
 //  - 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:
+// 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(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
 {
@@ -57,45 +87,155 @@ 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>;
 
-    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];
+    }
+    [[nodiscard]] HitT const& last() const
+    {
+        TG_ASSERT(_size > 0);
+        return _hit[_size - 1];
+    }
 
-    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)
+    hits() = default;
+    hits(HitT* hits, int size) : _size(size)
     {
         for (auto i = 0; i < size; ++i)
             _hit[i] = hits[i];
     }
+    template <typename... HitTs>
+    hits(HitTs... hits) : _size(sizeof...(HitTs)), _hit{hits...}
+    {
+    }
 
 private:
     int _size = 0;
-    HitT _hit[MaxHits];
+    HitT _hit[MaxHits] = {};
+};
+
+/// describes a continuous interval on a line or ray between start and end
+template <class ScalarT>
+struct hit_interval
+{
+    static constexpr bool is_hit_interval = true; // tag
+
+    ScalarT start;
+    ScalarT end;
+
+    [[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 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};
+}
+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]...);
+}
+
+// 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) || ...);
+}
+
+// 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&>()));
+}
+
+
 // ====================================== 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)
+    -> 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 <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) -> decltype(closest_intersection_parameter(a, b).has_value())
+{
+    return closest_intersection_parameter(a, b).has_value();
+}
+
+// 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();
 }
@@ -132,29 +272,81 @@ template <class A, class B>
     return {};
 }
 
-// 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>>
+// 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_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 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>>
+// 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_hit_interval, // condition
+    optional<segment<object_traits<A>::domain_dimension, typename object_traits<A>::scalar_t>>> // return type
+{
+    auto ts = intersection_parameter(a, b);
+    if (ts.has_value())
+        return {{a[ts.value().start], a[ts.value().end]}};
+    return {};
+}
+
+// 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>
 {
-    auto hits = intersection_parameter(r, obj);
+    const auto hits = intersection_parameter(a, b);
     if (hits.any())
-        return hits[0];
+        return hits.first();
     return {};
 }
 
+// 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)>
+{
+    const auto hits = intersection_parameter(a, b);
+    if (hits.has_value())
+        return hits.value().start;
+    return {};
+}
+
+// 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<hit_interval<ScalarT>>>
+{
+    const hits<2, ScalarT> inter = intersection_parameter(l, boundary_of(obj));
+
+    if (inter.size() == 2)
+        return {{inter[0], inter[1]}};
+
+    if constexpr (object_traits<Obj>::is_finite)
+    {
+        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() == 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
 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)
@@ -174,218 +366,561 @@ constexpr optional<pos<D, ScalarT>> intersection(Obj const& obj, pos<D, ScalarT>
 }
 
 
-// ====================================== Ray - Object Intersections ======================================
+// ====================================== Ray Intersections from Line Intersections ======================================
 
-// ray - plane
-template <int D, class ScalarT>
-[[nodiscard]] constexpr optional<ScalarT> intersection_parameter(ray<D, ScalarT> const& r, plane<D, ScalarT> const& p)
+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_hits, intersection_parameter(inf_of(ray), obj))
 {
-    // if plane normal and ray direction are parallel there is no intersection
-    auto dotND = dot(p.normal, r.dir);
-    if (dotND == ScalarT(0))
+    const auto inter = intersection_parameter(inf_of(ray), obj);
+    constexpr auto maxHits = inter.max_hits;
+    if (!inter.any() || inter.last() < ScalarT(0))
         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>
+    if constexpr (maxHits == 2)
+    {
+        if (inter.first() < ScalarT(0))
+            return {inter[1]};
+    }
+    else
+        static_assert(maxHits == 1, "Only up to two intersections implemented");
 
-    auto t = (p.dis - dot(p.normal, vec<D, ScalarT>(r.origin))) / dotND;
+    return inter;
+}
+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_hit_interval, intersection_parameter(inf_of(ray), obj))
+{
+    const auto inter = intersection_parameter(inf_of(ray), obj);
 
-    // check whether plane lies behind ray
-    if (t < 0)
+    if (!inter.has_value())
+        return inter;
+
+    auto interval = inter.value();
+    if (interval.end < ScalarT(0))
         return {};
 
-    return t;
+    TG_ASSERT((interval.start <= ScalarT(0)) == contains(obj, ray.origin));
+
+    interval.start = max(interval.start, ScalarT(0));
+    return interval;
 }
 
-// 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)
+// ====================================== Line - Object Intersections ======================================
+
+// line - point
+template <class ScalarT>
+[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<1, ScalarT> const& l, pos<1, ScalarT> const& p)
 {
-    auto t = dot(s.center - r.origin, r.dir);
+    return coordinates(l, p);
+}
 
-    auto d_sqr = distance_sqr(r[t], s.center);
-    auto r_sqr = s.radius * s.radius;
-    if (d_sqr > r_sqr)
+// line - line
+template <class ScalarT>
+[[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);
+    auto t = inverse(M) * (l1.pos - l0.pos);
+    return t.x;
+}
+
+// line - ray
+template <class ScalarT>
+[[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);
+    if (t.y < ScalarT(0))
         return {};
+    return t.x;
+}
 
-    auto dt = sqrt(r_sqr - d_sqr);
+// line - segment
+template <class ScalarT>
+[[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);
+    if (t.y < ScalarT(0) || t.y > ScalarT(1))
+        return {};
+    return t.x;
+}
 
-    if (t - dt >= 0)
-        return t - dt;
-    if (t + dt >= 0)
-        return t + dt;
+// line - plane
+template <int D, class ScalarT>
+[[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
+        return {};
 
-    return {};
+    // <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;
 }
 
-// returns intersection point(s) of ray and sphere
+// line - halfspace
 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 optional<hit_interval<ScalarT>> intersection_parameter(line<D, ScalarT> const& l, halfspace<D, ScalarT> const& h)
 {
-    auto t = dot(s.center - r.origin, r.dir);
+    const auto dotND = dot(h.normal, l.dir);
+    const auto dist = signed_distance(l.pos, h);
 
-    auto d_sqr = distance_sqr(r[t], s.center);
-    auto r_sqr = s.radius * s.radius;
-    if (d_sqr > r_sqr)
-        return {};
+    if (dotND == ScalarT(0)) // if halfspace normal and line direction are orthogonal, there is no intersection
+    {
+        if (dist <= ScalarT(0))
+            return {{tg::min<ScalarT>(), tg::max<ScalarT>()}}; // completely contained
+        return {}; // completely outside
+    }
 
-    auto dt = sqrt(r_sqr - d_sqr);
+    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)
+{
+    // 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
+    const auto dotND = dot(h.normal, r.dir);
+    if (dotND >= ScalarT(0))
+        return {};
 
-    ScalarT hits[2];
+    return -dist / dotND;
+}
 
-    if (t - dt >= 0)
+// line - aabb
+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)
     {
-        hits[0] = t - dt;
-        hits[1] = t + dt;
-        return {hits, 2};
+        if (abs(l.dir[i]) > ScalarT(100) * tg::epsilon<ScalarT>)
+        {
+            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 (l.pos[i] < b.min[i] || l.pos[i] > b.max[i])
+            return {}; // ray parallel to this axis and outside of the aabb
     }
 
-    if (t + dt >= 0)
+    // no intersection if line misses the aabb
+    if (tFirst > tSecond)
+        return {};
+
+    return {tFirst, tSecond};
+}
+
+// line - box
+template <int D, class ScalarT>
+[[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;
+    auto tFirst = tg::min<ScalarT>();
+    auto tSecond = tg::max<ScalarT>();
+    for (auto i = 0; i < D; ++i)
     {
-        hits[0] = t + dt;
-        return {hits, 1};
+        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;
+            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
     }
 
-    return {};
-}
+    // no intersection if line misses the box
+    if (tFirst > tSecond)
+        return {};
 
-// ray - tube
+    return {tFirst, tSecond};
+}
 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 hits<1, ScalarT> intersection_parameter(line<3, ScalarT> const& l, box<2, ScalarT, 3> const& b)
 {
-    auto cdir = direction(c);
-    auto cosA = dot(cdir, r.dir);
-    auto sinA_sqr = 1 - cosA * cosA;
+    const auto t = intersection_parameter(l, plane<3, ScalarT>(normal_of(b), b.center));
+    if (!t.any()) // line parallel to box plane
+        return {};
 
-    if (sinA_sqr <= 0)
+    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 {};
 
-    // 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;
+    return t;
+}
 
-    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;
+// line - disk
+template <class ScalarT>
+[[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
+        return {};
 
-    if (line_ray_dist_sqr > cyl_radius_sqr)
+    const auto p = l[t.first()];
+    if (distance_sqr(p, d.center) > d.radius * d.radius)
         return {};
 
-    // radius in 2D slice
-    auto r_2D = sqrt(cyl_radius_sqr - line_ray_dist_sqr);
+    return t;
+}
 
-    // 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;
+// line - sphere
+template <int D, class ScalarT>
+[[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);
 
-    // 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);
+    auto d_sqr = distance_sqr(l[t], s.center);
+    auto r_sqr = s.radius * s.radius;
+    if (d_sqr > r_sqr)
+        return {};
 
-    ScalarT hits[2];
-    int hit_cnt = 0;
+    auto dt = sqrt(r_sqr - d_sqr);
+    return {t - dt, t + dt};
+}
 
-    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;
+// line - hemisphere
+template <int D, class ScalarT>
+[[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 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(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, l[hit]))
+            hits[numHits++] = hit;
+
+    return {hits, numHits};
+}
+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(lr, h).any();
+    else
+        return detail::intersects_any(lr, caps_of(h), boundary_no_caps_of(h));
+}
 
-    return {hits, hit_cnt};
+// line - capsule
+template <class ScalarT>
+[[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);
+    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 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(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 - disk
+// line - cylinder
 template <class ScalarT>
-[[nodiscard]] constexpr optional<ScalarT> intersection_parameter(ray<3, ScalarT> const& r, sphere<2, ScalarT, 3> const& d)
+[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, cylinder_boundary<3, ScalarT> const& c)
 {
-    auto const t = intersection_parameter(r, plane<3, ScalarT>(d.normal, d.center));
-    if (!t.has_value())
-        return {};
+    const auto caps = caps_of(c);
+    return detail::merge_hits(l, caps[0], caps[1], boundary_no_caps_of(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(lr, c).any();
+    else
+    {
+        const auto caps = caps_of(c);
+        return detail::intersects_any(lr, caps[0], caps[1], boundary_no_caps_of(c));
+    }
+}
+template <class ScalarT>
+[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, cylinder_boundary_no_caps<3, ScalarT> const& c)
+{
+    auto const infInter = intersection_parameter(l, inf_of(c));
+    if (!infInter.any())
+        return infInter;
+
+    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);
+
+    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];
+
+    return {hits, numHits};
+}
 
-    auto const p = r[t.value()];
-    if (distance_sqr(p, d.center) > d.radius * d.radius)
+// 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;
+
+    if (sinA_sqr <= 0)
+        return {}; // line and cylinder are parallel
+
+    // 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;
+
+    auto const line_axis_dist_sqr = distance_sqr(l[tLine], c.axis[tAxis]);
+    auto const cyl_radius_sqr = c.radius * c.radius;
+
+    if (cyl_radius_sqr < line_axis_dist_sqr)
+        return {}; // line misses the 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)
+{
+    const auto n = perpendicular(c.axis.dir);
+    const auto d = dot(l.dir, n);
+    if (d == ScalarT(0)) // line parallel to inf_cylinder
         return {};
 
-    return t;
+    const auto dist = dot(c.axis.pos - l.pos, n);
+    const auto [tMin, tMax] = minmax((dist - c.radius) / d, (dist + c.radius) / d);
+    return {tMin, tMax};
 }
 
-// ray - quadric (as an isosurface, not error quadric)
+// 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 ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, quadric<3, ScalarT> const& q)
+[[nodiscard]] constexpr hits<2, ScalarT> intersection_parameter(line<3, ScalarT> const& l, inf_cone_boundary<3, ScalarT> const& ic)
 {
-    const auto Ad = q.A() * r.dir;
-    const auto p = r.origin;
+    // 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);
+    if (!inter.any())
+        return inter;
+
+    // exclude intersections with mirrored cone
+    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:
+    // 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, ic.opening_dir) >= ScalarT(0))
+        hits[numHits++] = inter[1];
+
+    return {hits, numHits};
+}
 
-    // 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;
+// 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};
+}
 
-    // 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
+// 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_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)
+{
+    return detail::merge_hits(l, py.base, boundary_no_caps_of(py));
+}
 
-    const auto sqrtD = sqrt(discriminant);
-    const auto t1 = (-b - sqrtD) / (ScalarT(2) * a);
-    const auto t2 = (-b + sqrtD) / (ScalarT(2) * a);
+// line - triangle2
+template <class ScalarT>
+[[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>();
+    auto numIntersections = 0;
+    for (const auto& edge : edges_of(t))
+    {
+        const auto inter = intersection_parameter(l, edge);
+        if (inter.any())
+        {
+            numIntersections++;
+            closestIntersection = min(closestIntersection, inter.first());
+            furtherIntersection = max(furtherIntersection, inter.first());
+        }
+    }
 
-    auto tMin = min(t1, t2);
-    auto tMax = max(t1, t2);
+    if (numIntersections == 0)
+        return {}; // No intersection
 
-    ScalarT hits[2];
+    TG_ASSERT(numIntersections == 2);
+    return {{closestIntersection, furtherIntersection}};
+}
 
-    if (tMin >= ScalarT(0))
-    {
-        hits[0] = tMin;
-        hits[1] = tMax;
-        return {hits, 2};
-    }
+// line - triangle3
+template <class ScalarT>
+[[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>)
+{
+    auto e1 = t.pos1 - t.pos0;
+    auto e2 = t.pos2 - t.pos0;
+
+    auto pvec = tg::cross(l.dir, e2);
+    auto det = dot(pvec, e1);
 
-    if (tMax >= ScalarT(0))
+    if (det < ScalarT(0))
     {
-        hits[0] = tMax;
-        return {hits, 1};
+        std::swap(e1, e2);
+        pvec = tg::cross(l.dir, e2);
+        det = -det;
     }
 
-    return {};
+    if (det < eps)
+        return {};
+
+    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(l.dir, qvec);
+    if (v < ScalarT(0) || v + u > det)
+        return {};
+
+    return (ScalarT(1) / det) * dot(e2, qvec);
 }
 
-// ray - cylinder
+// line - ellipse
 template <class ScalarT>
-[[nodiscard]] constexpr optional<ScalarT> closest_intersection_parameter(ray<3, ScalarT> const& r, cylinder<3, ScalarT> const& c)
+[[nodiscard]] constexpr hits<1, ScalarT> intersection_parameter(line<3, ScalarT> const& l, ellipse<2, ScalarT, 3> const& e)
 {
-    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));
-
-    optional<ScalarT> t;
-
-    if (t_cyl.has_value())
-        t = t_cyl.value();
+    const auto t = intersection_parameter(l, plane<3, ScalarT>(normal_of(e), e.center));
+    if (!t.any()) // line parallel to ellipse plane
+        return {};
 
-    if (t_cap0.has_value() && (!t.has_value() || t.value() > t_cap0.value()))
-        t = t_cap0;
+    // 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 (t_cap1.has_value() && (!t.has_value() || t.value() > t_cap1.value()))
-        t = t_cap1;
+    if (pow2(x / a) + pow2(y / b) <= ScalarT(1))
+        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;
+    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)
+    {
+        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));
+}
 
-    TG_INTERNAL_ASSERT(!t.has_value() || t.value() >= 0);
+// 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)
+{
+    const auto Ad = q.A() * l.dir;
+    const auto p = l.pos;
 
-    return t;
+    // 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(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;
+    return detail::solve_quadratic(a, b, c);
 }
 
 
+// ====================================== 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)
 template <class ScalarT>
@@ -548,14 +1083,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)
 {
@@ -588,12 +1115,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 - pos<3, ScalarT>::zero, 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)
 {
@@ -607,112 +1128,29 @@ 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;
 }
 
-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)
@@ -753,47 +1191,20 @@ 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 <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)
 {
+    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);
@@ -827,11 +1238,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>
@@ -840,7 +1246,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);
@@ -861,12 +1267,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
@@ -937,10 +1339,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
diff --git a/src/typed-geometry/functions/objects/normal.hh b/src/typed-geometry/functions/objects/normal.hh
index dd09e8bd01ae2ce448901fcd70cd13d2e578f1b3..9096e06b07bf5d19facfe87b1122ea21fb86cad6 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)
 {
diff --git a/src/typed-geometry/functions/objects/project.hh b/src/typed-geometry/functions/objects/project.hh
index 1794cbe9ffa49fccd18fe8b8dc7a90cc8dba3d72..26747e81aef4af89f2ad9621f71f3a55eca0eb39 100644
--- a/src/typed-geometry/functions/objects/project.hh
+++ b/src/typed-geometry/functions/objects/project.hh
@@ -30,12 +30,6 @@
 
 namespace tg
 {
-template <int D, class ScalarT>
-[[nodiscard]] constexpr pos<D, ScalarT> project(pos<D, ScalarT> const&, pos<D, ScalarT> const& target_p)
-{
-    return target_p;
-}
-
 // ============== project to plane ==============
 
 template <int D, class ScalarT>
@@ -66,6 +60,15 @@ template <int D, class ScalarT>
 }
 
 
+// ============== project to point ==============
+
+template <int D, class ScalarT>
+[[nodiscard]] constexpr pos<D, ScalarT> project([[maybe_unused]] pos<D, ScalarT> const& p, pos<D, ScalarT> const& q)
+{
+    return q;
+}
+
+
 // ============== project to line / ray / segment ==============
 
 template <int D, class ScalarT>
diff --git a/src/typed-geometry/types/objects/aabb.hh b/src/typed-geometry/types/objects/aabb.hh
index b0e9121bc980d345c6c27313a7ba9cf758f6b37c..5afc17dccd790b5addf61ecc0298aa457cbe7c3d 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>;
 
diff --git a/src/typed-geometry/types/objects/inf_cone.hh b/src/typed-geometry/types/objects/inf_cone.hh
index 6e171639483955128def72e8e9fa8217b3fa968f..1af69f0c3f8f9e71b26e1ce6a52e4c4e3bca592a 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 8d1aa96307109ec7a913d3966a2ff73797a1b91b..259002228b95a28219c545a70f5b5417950f763e 100644
--- a/src/typed-geometry/types/objects/inf_cylinder.hh
+++ b/src/typed-geometry/types/objects/inf_cylinder.hh
@@ -8,16 +8,22 @@
 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_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 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>;
@@ -29,6 +35,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>;
diff --git a/src/typed-geometry/types/objects/line.hh b/src/typed-geometry/types/objects/line.hh
index d0f84b5724a9d0539c89f358dbb2b95f110bca6a..80071c7ea1e8a39912d367aced55a0cf767692d8 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,7 +45,8 @@ using uline4 = line<4, u32>;
 template <int D, class ScalarT>
 struct line
 {
-    using dir_t = tg::dir<D, ScalarT>;
+    using scalar_t = 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;
diff --git a/src/typed-geometry/types/objects/ray.hh b/src/typed-geometry/types/objects/ray.hh
index 402f3d004383884da68fbe33b03852c887a92fd8..f0b85dfe5bd76b5abf15127b1be9cc9c9332e2ce 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,7 +43,8 @@ using uray4 = ray<4, u32>;
 template <int D, class ScalarT>
 struct ray
 {
-    using dir_t = tg::dir<D, ScalarT>;
+    using scalar_t = 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;
diff --git a/src/typed-geometry/types/objects/segment.hh b/src/typed-geometry/types/objects/segment.hh
index 8f6f91a0e731c68d01583fcedc04a2143c87703e..d80e584403a47e153d7aa95884ced2eb24e8d476 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;
diff --git a/src/typed-geometry/types/objects/tube.hh b/src/typed-geometry/types/objects/tube.hh
index bc90acd801ca53511a615a4127063b042cfd110d..16e1f543933196fd66f982694b281b77f5bd9299 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>;
 }