WIP: moving to new traits arch

src/typed-geometry/detail/constructors-fun.hh

@@ -8,11 +8,11 @@
 #include <typed-geometry/types/vec.hh>
 
 #include <typed-geometry/types/objects/box.hh>
-#include <typed-geometry/types/objects/halfspace.hh>
-#include <typed-geometry/types/objects/plane.hh>
-#include <typed-geometry/types/objects/line.hh>
 #include <typed-geometry/types/objects/cone.hh>
+#include <typed-geometry/types/objects/halfspace.hh>
 #include <typed-geometry/types/objects/inf_cone.hh>
+#include <typed-geometry/types/objects/line.hh>
+#include <typed-geometry/types/objects/plane.hh>
 
 #include <typed-geometry/detail/operators/ops_pos.hh>
 #include <typed-geometry/functions/vector/dot.hh>
@@ -21,8 +21,8 @@
 
 namespace tg
 {
-template <int D, class ScalarT>
-constexpr box<D, ScalarT>::box(aabb<D, ScalarT> const& b)
+template <int D, class ScalarT, class TraitsT>
+constexpr box<D, ScalarT, D, TraitsT>::box(aabb<D, ScalarT, TraitsT> const& b)
 {
     static_assert(is_floating_point<ScalarT>, "cannot be guaranteed for integers");
     auto half_e = (b.max - b.min) / ScalarT(2);
@@ -50,8 +50,8 @@ constexpr line<D, ScalarT> line<D, ScalarT>::from_points(pos_t a, pos_t b)
     return line(a, normalize(b - a));
 }
 
-template <int D, class ScalarT>
-constexpr inf_cone<D, ScalarT>::inf_cone(cone<D, ScalarT> c)
+template <int D, class ScalarT, class TraitsT>
+constexpr inf_cone<D, ScalarT, TraitsT>::inf_cone(cone<D, ScalarT, TraitsT> const& c)
 {
     apex = c.base.center + c.height * c.base.normal;
     opening_dir = -c.base.normal;

src/typed-geometry/detail/operators/subscript.hh

@@ -33,8 +33,8 @@ template <int D, class ScalarT>
     return this->pos + this->dir * t;
 }
 
-template <int D, class ScalarT>
-[[nodiscard]] constexpr pos<D, ScalarT> aabb<D, ScalarT>::operator[](comp<D, ScalarT> const& c) const
+template <int D, class ScalarT, class TraitsT>
+[[nodiscard]] constexpr pos<D, ScalarT> aabb<D, ScalarT, TraitsT>::operator[](comp<D, ScalarT> const& c) const
 {
     return this->min + (this->max - this->min) * size(c);
 }

src/typed-geometry/detail/special_values.hh

@@ -321,21 +321,22 @@ constexpr mat<C, R, ScalarT> mat<C, R, ScalarT>::diag(vec<detail::min(C, R), Sca
 // Objects
 
 
-template <int D, class ScalarT>
-const aabb<D, ScalarT> aabb<D, ScalarT>::minus_one_to_one = {tg::pos<D, ScalarT>(ScalarT(-1)), tg::pos<D, ScalarT>(ScalarT(1))};
-template <int D, class ScalarT>
-const aabb<D, ScalarT> aabb<D, ScalarT>::unit_from_zero = {tg::pos<D, ScalarT>(ScalarT(0)), tg::pos<D, ScalarT>(ScalarT(1))};
-template <int D, class ScalarT>
-const aabb<D, ScalarT> aabb<D, ScalarT>::unit_centered = {tg::pos<D, ScalarT>(ScalarT(-0.5)), tg::pos<D, ScalarT>(ScalarT(0.5))};
-
-template <int D, class ScalarT>
-const box<D, ScalarT> box<D, ScalarT>::minus_one_to_one = {tg::pos<D, ScalarT>::zero, tg::mat<D, D, ScalarT>::identity};
-template <int D, class ScalarT>
-const box<D, ScalarT> box<D, ScalarT>::unit_from_zero = {tg::pos<D, ScalarT>(ScalarT(0.5)), tg::mat<D, D, ScalarT>::diag(ScalarT(0.5))};
-template <int D, class ScalarT>
-const box<D, ScalarT> box<D, ScalarT>::unit_centered = {tg::pos<D, ScalarT>::zero, tg::mat<D, D, ScalarT>::diag(ScalarT(0.5))};
-
-template <int D, class ScalarT>
-const sphere<D, ScalarT> sphere<D, ScalarT>::unit = {pos<D, ScalarT>::zero, ScalarT(1)};
+template <int D, class ScalarT, class TraitsT>
+const aabb<D, ScalarT, TraitsT> aabb<D, ScalarT, TraitsT>::minus_one_to_one = {tg::pos<D, ScalarT>(ScalarT(-1)), tg::pos<D, ScalarT>(ScalarT(1))};
+template <int D, class ScalarT, class TraitsT>
+const aabb<D, ScalarT, TraitsT> aabb<D, ScalarT, TraitsT>::unit_from_zero = {tg::pos<D, ScalarT>(ScalarT(0)), tg::pos<D, ScalarT>(ScalarT(1))};
+template <int D, class ScalarT, class TraitsT>
+const aabb<D, ScalarT, TraitsT> aabb<D, ScalarT, TraitsT>::unit_centered = {tg::pos<D, ScalarT>(ScalarT(-0.5)), tg::pos<D, ScalarT>(ScalarT(0.5))};
+
+template <int D, class ScalarT, class TraitsT>
+const box<D, ScalarT, D, TraitsT> box<D, ScalarT, D, TraitsT>::minus_one_to_one = {tg::pos<D, ScalarT>::zero, tg::mat<D, D, ScalarT>::identity};
+template <int D, class ScalarT, class TraitsT>
+const box<D, ScalarT, D, TraitsT> box<D, ScalarT, D, TraitsT>::unit_from_zero
+    = {tg::pos<D, ScalarT>(ScalarT(0.5)), tg::mat<D, D, ScalarT>::diag(ScalarT(0.5))};
+template <int D, class ScalarT, class TraitsT>
+const box<D, ScalarT, D, TraitsT> box<D, ScalarT, D, TraitsT>::unit_centered = {tg::pos<D, ScalarT>::zero, tg::mat<D, D, ScalarT>::diag(ScalarT(0.5))};
+
+template <int D, class ScalarT, class TraitsT>
+const sphere<D, ScalarT, D, TraitsT> sphere<D, ScalarT, D, TraitsT>::unit = {pos<D, ScalarT>::zero, ScalarT(1)};
 
 } // namespace tg

src/typed-geometry/feature/objects.hh

@@ -4,6 +4,7 @@
 
 #include <typed-geometry/functions/objects/aabb.hh>
 #include <typed-geometry/functions/objects/area.hh>
+#include <typed-geometry/functions/objects/boundary.hh>
 #include <typed-geometry/functions/objects/bounds.hh>
 #include <typed-geometry/functions/objects/centroid.hh>
 #include <typed-geometry/functions/objects/closest_points.hh>

src/typed-geometry/functions/objects/boundary.hh

+#pragma once
+
+#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/hemisphere.hh>
+#include <typed-geometry/types/objects/inf_cone.hh>
+#include <typed-geometry/types/objects/inf_cylinder.hh>
+#include <typed-geometry/types/objects/pyramid.hh>
+#include <typed-geometry/types/objects/sphere.hh>
+
+// tg::boundary_of(obj) converts an object to its boundary (identity for objects that are already boundaries)
+
+namespace tg
+{
+template <int D, class ScalarT, class TraitsT>
+[[nodiscard]] constexpr aabb_boundary<D, ScalarT> boundary_of(aabb<D, ScalarT, TraitsT> const& v)
+{
+    return {v.min, v.max};
+}
+template <int D, class ScalarT, class TraitsT>
+[[nodiscard]] constexpr box_boundary<D, ScalarT> boundary_of(box<D, ScalarT, D, TraitsT> const& v)
+{
+    return {v.center, v.half_extents};
+}
+template <class ScalarT, class TraitsT>
+[[nodiscard]] constexpr box_boundary<2, ScalarT, 3> boundary_of(box<2, ScalarT, 3, TraitsT> const& v)
+{
+    return {v.center, v.half_extents, v.normal};
+}
+template <int D, class ScalarT, class TraitsT>
+[[nodiscard]] constexpr hemisphere_boundary<D, ScalarT> boundary_of(hemisphere<D, ScalarT, TraitsT> const& v)
+{
+    return {v.center, v.radius, v.normal};
+}
+template <int D, class ScalarT, class TraitsT>
+[[nodiscard]] constexpr capsule_boundary<D, ScalarT> boundary_of(capsule<D, ScalarT, TraitsT> const& v)
+{
+    return {v.axis, v.radius};
+}
+template <int D, class ScalarT, class TraitsT>
+[[nodiscard]] constexpr cylinder_boundary<D, ScalarT> boundary_of(cylinder<D, ScalarT, TraitsT> const& v)
+{
+    return {v.axis, v.radius};
+}
+template <int D, class ScalarT, class TraitsT>
+[[nodiscard]] constexpr inf_cone_boundary<D, ScalarT> boundary_of(inf_cone<D, ScalarT, TraitsT> const& v)
+{
+    return {v.apex, v.opening_dir, v.opening_angle};
+}
+template <int D, class ScalarT, class TraitsT>
+[[nodiscard]] constexpr inf_cylinder_boundary<D, ScalarT> boundary_of(inf_cylinder<D, ScalarT, TraitsT> const& v)
+{
+    return {v.axis, v.radius};
+}
+template <class BaseT, class TraitsT>
+[[nodiscard]] constexpr pyramid_boundary<BaseT> boundary_of(pyramid<BaseT, TraitsT> const& v)
+{
+    return {v.base, v.height};
+}
+template <int D, class ScalarT, class TraitsT>
+[[nodiscard]] constexpr sphere_boundary<D, ScalarT> boundary_of(sphere<D, ScalarT, D, TraitsT> const& v)
+{
+    return {v.center, v.radius};
+}
+template <class ScalarT, class TraitsT>
+[[nodiscard]] constexpr sphere_boundary<2, ScalarT, 3> boundary_of(sphere<2, ScalarT, 3, TraitsT> const& v)
+{
+    return {v.center, v.radius, v.normal};
+}
+
+// === no caps versions ===
+
+template <int D, class ScalarT>
+[[nodiscard]] constexpr cylinder_boundary_no_caps<D, ScalarT> boundary_of(cylinder_boundary_no_caps<D, ScalarT> v)
+{
+    return v;
+}
+template <int D, class ScalarT>
+[[nodiscard]] constexpr hemisphere_boundary_no_caps<D, ScalarT> boundary_of(hemisphere_boundary_no_caps<D, ScalarT> v)
+{
+    return v;
+}
+template <class BaseT>
+[[nodiscard]] constexpr pyramid_boundary_no_caps<BaseT> boundary_of(pyramid_boundary_no_caps<BaseT> v)
+{
+    return v;
+}
+
+}

src/typed-geometry/functions/objects/contains.hh

@@ -3,7 +3,6 @@
 #include <cmath>
 
 #include <typed-geometry/types/objects/aabb.hh>
-#include <typed-geometry/types/objects/ball.hh>
 #include <typed-geometry/types/objects/box.hh>
 #include <typed-geometry/types/objects/cylinder.hh>
 #include <typed-geometry/types/objects/sphere.hh>
@@ -29,13 +28,13 @@ namespace tg
 {
 // default implementation if distance(a, pos) is available
 template <class A, int D, class ScalarT>
-[[nodiscard]] constexpr auto contains(A const& a, pos<D, ScalarT> const& p, ScalarT eps = ScalarT(0)) -> decltype(ScalarT(distance(a, p)), false)
+[[nodiscard]] constexpr auto contains(A const& a, pos<D, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0)) -> decltype(ScalarT(distance(a, p)), false)
 {
     return distance(a, p) <= eps;
 }
 
 template <int D, class ScalarT>
-[[nodiscard]] constexpr bool contains(pos<D, ScalarT> const& b, pos<D, ScalarT> const& o, ScalarT eps = ScalarT(0))
+[[nodiscard]] constexpr bool contains(pos<D, ScalarT> const& b, pos<D, ScalarT> const& o, dont_deduce<ScalarT> eps = ScalarT(0))
 {
     if (eps > ScalarT(0))
         return distance_sqr(b, o) < eps * eps;
@@ -43,30 +42,30 @@ template <int D, class ScalarT>
 }
 
 template <class ScalarT>
-[[nodiscard]] constexpr bool contains(aabb<1, ScalarT> const& b, ScalarT const& o, ScalarT eps = ScalarT(0))
+[[nodiscard]] constexpr bool contains(aabb<1, ScalarT> const& b, ScalarT const& o, dont_deduce<ScalarT> eps = ScalarT(0))
 {
     return b.min.x - eps <= o && o <= b.max.x + eps;
 }
 template <class ScalarT>
-[[nodiscard]] constexpr bool contains(aabb<1, ScalarT> const& b, pos<1, ScalarT> const& o, ScalarT eps = ScalarT(0))
+[[nodiscard]] constexpr bool contains(aabb<1, ScalarT> const& b, pos<1, ScalarT> const& o, dont_deduce<ScalarT> eps = ScalarT(0))
 {
     return b.min.x - eps <= o.x && o.x <= b.max.x + eps;
 }
 template <class ScalarT>
-[[nodiscard]] constexpr bool contains(aabb<2, ScalarT> const& b, pos<2, ScalarT> const& o, ScalarT eps = ScalarT(0))
+[[nodiscard]] constexpr bool contains(aabb<2, ScalarT> const& b, pos<2, ScalarT> const& o, dont_deduce<ScalarT> eps = ScalarT(0))
 {
     return b.min.x - eps <= o.x && o.x <= b.max.x + eps && //
            b.min.y - eps <= o.y && o.y <= b.max.y + eps;
 }
 template <class ScalarT>
-[[nodiscard]] constexpr bool contains(aabb<3, ScalarT> const& b, pos<3, ScalarT> const& o, ScalarT eps = ScalarT(0))
+[[nodiscard]] constexpr bool contains(aabb<3, ScalarT> const& b, pos<3, ScalarT> const& o, dont_deduce<ScalarT> eps = ScalarT(0))
 {
     return b.min.x - eps <= o.x && o.x <= b.max.x + eps && //
            b.min.y - eps <= o.y && o.y <= b.max.y + eps && //
            b.min.z - eps <= o.z && o.z <= b.max.z + eps;
 }
 template <class ScalarT>
-[[nodiscard]] constexpr bool contains(aabb<4, ScalarT> const& b, pos<4, ScalarT> const& o, ScalarT eps = ScalarT(0))
+[[nodiscard]] constexpr bool contains(aabb<4, ScalarT> const& b, pos<4, ScalarT> const& o, dont_deduce<ScalarT> eps = ScalarT(0))
 {
     return b.min.x - eps <= o.x && o.x <= b.max.x + eps && //
            b.min.y - eps <= o.y && o.y <= b.max.y + eps && //
@@ -75,7 +74,7 @@ template <class ScalarT>
 }
 
 template <int D, class ScalarT>
-[[nodiscard]] constexpr bool contains(box<D, ScalarT> const& b, pos<D, ScalarT> const& o, ScalarT eps = ScalarT(0))
+[[nodiscard]] constexpr bool contains(box<D, ScalarT> const& b, pos<D, ScalarT> const& o, dont_deduce<ScalarT> eps = ScalarT(0))
 {
     auto r = o - b.center;
     // TODO: unroll
@@ -86,30 +85,21 @@ template <int D, class ScalarT>
 }
 
 template <int D, class ScalarT>
-[[nodiscard]] constexpr bool contains(aabb<D, ScalarT> const& b, aabb<D, ScalarT> const& o, ScalarT eps = ScalarT(0))
+[[nodiscard]] constexpr bool contains(aabb<D, ScalarT> const& b, aabb<D, ScalarT> const& o, dont_deduce<ScalarT> eps = ScalarT(0))
 {
     return contains(b, o.min, eps) && contains(b, o.max, eps);
 }
 
 template <int D, class ScalarT>
-[[nodiscard]] constexpr bool contains(ball<D, ScalarT> const& s, pos<D, ScalarT> const& p, ScalarT eps = ScalarT(0))
+[[nodiscard]] constexpr bool contains(sphere<D, ScalarT> const& s, pos<D, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0))
 {
     auto r = s.radius + eps;
     return distance_sqr(s.center, p) <= r * r;
 }
 
-template <int D, class ScalarT>
-[[nodiscard]] constexpr bool contains(sphere<D, ScalarT> const& s, pos<D, ScalarT> const& p, ScalarT eps = ScalarT(0))
-{
-    auto const d2 = distance_sqr(s.center, p);
-    auto const s2 = s.radius * s.radius;
-    auto const e2 = eps * eps;
-    return d2 - e2 <= s2 && s2 <= d2 + e2;
-}
-
 // Note that eps is used to compare 2D areas, not 1D lengths
 template <class ScalarT>
-[[nodiscard]] constexpr bool contains(triangle<2, ScalarT> const& t, pos<2, ScalarT> const& p, ScalarT eps = ScalarT(0))
+[[nodiscard]] constexpr bool contains(triangle<2, ScalarT> const& t, pos<2, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0))
 {
     auto pv0 = t.pos0 - p;
     auto pv1 = t.pos1 - p;
@@ -177,19 +167,7 @@ template <class ScalarT>
 }
 
 template <class ScalarT>
-[[nodiscard]] constexpr bool contains(circle<2, ScalarT> const& c, pos<2, ScalarT> const& p, ScalarT eps = ScalarT(0))
-{
-    return abs(distance_sqr(c.center, p) - c.radius * c.radius) <= eps * eps;
-}
-template <class ScalarT>
-[[nodiscard]] constexpr bool contains(disk<2, ScalarT> const& d, pos<2, ScalarT> const& p, ScalarT eps = ScalarT(0))
-{
-    auto r = d.radius + eps;
-    return distance_sqr(d.center, p) <= r * r;
-}
-
-template <class ScalarT>
-[[nodiscard]] constexpr bool contains(cone<3, ScalarT> const& c, pos<3, ScalarT> const& p, ScalarT eps = ScalarT(0))
+[[nodiscard]] constexpr bool contains(cone<3, ScalarT> const& c, pos<3, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0))
 {
     auto center = c.base.center - eps * c.base.normal;
 
@@ -204,7 +182,7 @@ template <class ScalarT>
 }
 
 template <class ScalarT>
-[[nodiscard]] constexpr bool contains(inf_cone<3, ScalarT> const& c, pos<3, ScalarT> const& p, ScalarT eps = ScalarT(0))
+[[nodiscard]] constexpr bool contains(inf_cone<3, ScalarT> const& c, pos<3, ScalarT> const& p, dont_deduce<ScalarT> eps = ScalarT(0))
 {
     auto apex = c.apex - (c.opening_dir * eps); // Shift apex outwards to add eps
     return angle_between(p - apex, c.opening_dir) <= c.opening_angle;

src/typed-geometry/functions/objects/direction.hh

@@ -5,7 +5,6 @@
 #include <typed-geometry/types/objects/line.hh>
 #include <typed-geometry/types/objects/ray.hh>
 #include <typed-geometry/types/objects/segment.hh>
-#include <typed-geometry/types/objects/tube.hh>
 
 #include <typed-geometry/functions/vector/normalize.hh>
 
@@ -40,11 +39,6 @@ template <int D, class ScalarT>
     return normalize(s.pos1 - s.pos0);
 }
 template <int D, class ScalarT>
-[[nodiscard]] constexpr dir<D, ScalarT> direction(tube<D, ScalarT> const& c)
-{
-    return direction(c.axis);
-}
-template <int D, class ScalarT>
 [[nodiscard]] constexpr dir<D, ScalarT> direction(cylinder<D, ScalarT> const& c)
 {
     return direction(c.axis);

src/typed-geometry/functions/objects/intersection.hh

 #pragma once
 
+#include <typed-geometry/detail/optional.hh>
 #include <typed-geometry/feature/assert.hh>
 #include <typed-geometry/functions/basic/scalar_math.hh>
-#include <typed-geometry/detail/optional.hh>
 
 #include <typed-geometry/types/objects/capsule.hh>
-#include <typed-geometry/types/objects/circle.hh>
 #include <typed-geometry/types/objects/cylinder.hh>
 #include <typed-geometry/types/objects/plane.hh>
 #include <typed-geometry/types/objects/ray.hh>
 #include <typed-geometry/types/objects/segment.hh>
 #include <typed-geometry/types/objects/sphere.hh>
 #include <typed-geometry/types/objects/triangle.hh>
-#include <typed-geometry/types/objects/tube.hh>
 
 #include <typed-geometry/functions/vector/cross.hh>
 #include <typed-geometry/functions/vector/dot.hh>
 #include <typed-geometry/functions/vector/length.hh>
 
-#include "closest_points.hh"
 #include "aabb.hh"
+#include "closest_points.hh"
 #include "contains.hh"
 #include "direction.hh"
 #include "normal.hh"
@@ -255,7 +253,7 @@ template <int D, class ScalarT>
 
 // ray - tube
 template <class ScalarT>
-[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter(ray<3, ScalarT> const& r, tube<3, ScalarT> const& c)
+[[nodiscard]] constexpr ray_hits<2, ScalarT> intersection_parameter_no_caps(ray<3, ScalarT> const& r, cylinder<3, ScalarT> const& c)
 {
     auto cdir = direction(c);
     auto cosA = dot(cdir, r.dir);
@@ -305,7 +303,7 @@ template <class ScalarT>
 
 // ray - disk
 template <class ScalarT>
-[[nodiscard]] constexpr optional<ScalarT> intersection_parameter(ray<3, ScalarT> const& r, disk<3, ScalarT> const& d)
+[[nodiscard]] constexpr optional<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())
@@ -366,9 +364,9 @@ 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, tube<3, ScalarT>(c.axis, c.radius));
-    auto const t_cap0 = intersection_parameter(r, disk<3, ScalarT>(c.axis.pos0, c.radius, dir));
-    auto const t_cap1 = intersection_parameter(r, disk<3, ScalarT>(c.axis.pos1, c.radius, dir));
+    auto const t_cyl = closest_intersection_parameter_no_caps(r, cylinder<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;
 
@@ -390,7 +388,7 @@ template <class ScalarT>
 // 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>
-[[nodiscard]] constexpr optional<circle<3, ScalarT>> intersection(sphere<3, ScalarT> const& a, sphere<3, ScalarT> const& b)
+[[nodiscard]] constexpr optional<sphere_boundary<2, ScalarT, 3>> intersection(sphere_boundary<3, ScalarT> const& a, sphere_boundary<3, ScalarT> const& b)
 {
     auto d2 = distance_sqr(a.center, b.center);
 
@@ -437,13 +435,13 @@ template <class ScalarT>
     auto irad = sqrt(ar2 - t * t * d2);
 
     // non-empty intersection (circle)
-    return tg::circle3{ipos, irad, dir<3, ScalarT>((b.center - a.center) / d)};
+    return sphere_boundary<2, ScalarT, 3>{ipos, irad, dir<3, ScalarT>((b.center - a.center) / d)};
 }
 
 // returns intersection points of two circles in 2D
 // for now does not work if circles are identical (result would be a circle2 again)
 template <class ScalarT>
-[[nodiscard]] constexpr optional<pair<pos<2, ScalarT>, pos<2, ScalarT>>> intersection(circle<2, ScalarT> const& a, circle<2, ScalarT> const& b)
+[[nodiscard]] constexpr optional<array<pos<2, ScalarT>, 2>> intersection(sphere_boundary<2, ScalarT> const& a, sphere_boundary<2, ScalarT> const& b)
 {
     if (a.center == b.center && a.radius == b.radius)
         return {}; // degenerate case
@@ -476,7 +474,7 @@ template <class ScalarT>
     auto p_above = p_between + h_by_d * a_to_b_swap;
     auto p_below = p_between - h_by_d * a_to_b_swap;
 
-    return tg::pair{p_above, p_below};
+    return array<pos<2, ScalarT>, 2>{p_above, p_below};
 }
 
 

src/typed-geometry/functions/objects/project.hh

@@ -14,6 +14,7 @@
 #include <typed-geometry/types/vec.hh>
 
 #include <typed-geometry/functions/matrix/inverse.hh>
+#include <typed-geometry/functions/vector/project.hh>
 
 #include "contains.hh"
 #include "coordinates.hh"
@@ -21,33 +22,37 @@
 
 namespace tg
 {
-// ============================== Vec Projections ==============================
+// ============== project to plane ==============
 
 template <int D, class ScalarT>
-[[nodiscard]] constexpr vec<D, ScalarT> project(vec<D, ScalarT> const& a, vec<D, ScalarT> const& b)
+[[nodiscard]] constexpr vec<D, ScalarT> project(vec<D, ScalarT> const& v, plane<D, ScalarT> const& pl)
 {
-    return b * dot(a, b) / dot(b, b);
+    return v - pl.normal * dot(v, pl.normal);
 }
 
 template <int D, class ScalarT>
-[[nodiscard]] constexpr vec<D, ScalarT> project(vec<D, ScalarT> const& a, dir<D, ScalarT> const& b)
+[[nodiscard]] constexpr vec<D, ScalarT> project(dir<D, ScalarT> const& v, plane<D, ScalarT> const& pl)
 {
-    return b * dot(a, b);
+    return v - pl.normal * dot(v, pl.normal);
 }
 
 template <int D, class ScalarT>
-[[nodiscard]] constexpr vec<D, ScalarT> project(vec<D, ScalarT> const& v, plane<D, ScalarT> const& pl)
+[[nodiscard]] constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, plane<D, ScalarT> const& pl)
 {
-    return v - pl.normal * dot(v, pl.normal);
+    return p - pl.normal * (dot(p, pl.normal) - pl.dis);
 }
 
+
+// ============== project to halfspace ==============
+
 template <int D, class ScalarT>
-[[nodiscard]] constexpr vec<D, ScalarT> project(dir<D, ScalarT> const& v, plane<D, ScalarT> const& pl)
+[[nodiscard]] constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, halfspace<D, ScalarT> const& pl)
 {
-    return v - pl.normal * dot(v, pl.normal);
+    return p - pl.normal * tg::max(ScalarT(0), dot(p, pl.normal) - pl.dis);
 }
 
-// ============================== Pos Projections ==============================
+
+// ============== project to line / ray / segment ==============
 
 template <int D, class ScalarT>
 [[nodiscard]] constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, line<D, ScalarT> const& l)
@@ -70,12 +75,18 @@ template <int D, class ScalarT>
     return s[t];
 }
 
+
+// ============== project to aabb ==============
+
 template <int D, class ScalarT>
 [[nodiscard]] constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, aabb<D, ScalarT> const& s)
 {
     return clamp(p, s.min, s.max);
 }
 
+
+// ============== project to box ==============
+
 template <int D, class ScalarT>
 [[nodiscard]] constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, box<D, ScalarT> const& b)
 {
@@ -83,17 +94,8 @@ template <int D, class ScalarT>
     return b.center + b.half_extents * vec<D, ScalarT>(project(pLocal, aabb<D, ScalarT>::minus_one_to_one));
 }
 
-template <int D, class ScalarT>
-[[nodiscard]] constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, plane<D, ScalarT> const& pl)
-{
-    return p - pl.normal * (dot(p, pl.normal) - pl.dis);
-}
 
-template <int D, class ScalarT>
-[[nodiscard]] constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, halfspace<D, ScalarT> const& pl)
-{
-    return p - pl.normal * tg::max(ScalarT(0), dot(p, pl.normal) - pl.dis);
-}
+// ============== project to triangle ==============
 
 template <class ScalarT>
 [[nodiscard]] constexpr pos<3, ScalarT> project(pos<3, ScalarT> const& p, triangle<3, ScalarT> const& t)
@@ -140,8 +142,20 @@ template <class ScalarT>
         return p2;
 }
 
+
+// ============== project to sphere ==============
+
+template <int D, class ScalarT>
+[[nodiscard]] constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, sphere<D, ScalarT> const& s)
+{
+    if (contains(s, p))
+        return p;
+
+    return project_to_boundary(p, s);
+}
+
 template <int D, class ScalarT>
-[[nodiscard]] constexpr pos<D, ScalarT> project(pos<D, ScalarT> const& p, sphere<D, ScalarT> const& sp)
+[[nodiscard]] constexpr pos<D, ScalarT> project_to_boundary(pos<D, ScalarT> const& p, sphere<D, ScalarT> const& sp)
 {
     auto dir_to_p = tg::normalize_safe(p - sp.center);
     if (is_zero_vector(dir_to_p))
@@ -149,6 +163,36 @@ template <int D, class ScalarT>
     return sp.center + dir_to_p * sp.radius;
 }
 
+template <class ScalarT>