Merge branch 'develop' of...

Merge branch 'develop' of https://www.graphics.rwth-aachen.de:9000/ptrettner/typed-geometry into develop

src/typed-geometry/detail/utility.hh

@@ -254,4 +254,16 @@ struct identity_fun
     }
 };
 
+namespace detail
+{
+// in detail ns to avoid ambiguity clash with std::swap
+template <class T>
+void swap(T& a, T& b)
+{
+    T tmp = static_cast<T&&>(a);
+    a = static_cast<T&&>(b);
+    b = static_cast<T&&>(tmp);
+}
+}
+
 } // namespace tg

src/typed-geometry/functions/angle.hh

@@ -24,14 +24,19 @@ TG_NODISCARD constexpr angle_t<fractional_result<ScalarT>> angle_between(vec<D,
 template <int D, class ScalarT>
 TG_NODISCARD constexpr angle_t<fractional_result<ScalarT>> angle_between(dir<D, ScalarT> const& a, dir<D, ScalarT> const& b)
 {
-    return acos(saturate(dot(a, b)));
+    constexpr auto lower = decltype (dot(normal(a), normal(b)))(-1);
+    constexpr auto upper = decltype (dot(normal(a), normal(b)))(1);
+    return acos(clamp(dot(a, b), lower, upper));
 }
 
 // returns the angle between any two objects with unambiguous normals. The result is in 0..pi (0°..180°)
 template <class A, class B>
 TG_NODISCARD constexpr auto angle_between(A const& a, B const& b) -> decltype(acos(dot(normal(a), normal(b))))
 {
-    return acos(saturate(dot(normal(a), normal(b))));
+    // TODO(ks): call to angle_between(dir, dir)?
+    constexpr auto lower = decltype (dot(normal(a), normal(b)))(-1);
+    constexpr auto upper = decltype (dot(normal(a), normal(b)))(1);
+    return acos(clamp(dot(normal(a), normal(b))), lower, upper);
 }
 
 // Returns the angle of a rotation of a towards b about the orthogonal_axis

src/typed-geometry/functions/intersection.hh

@@ -452,6 +452,30 @@ TG_NODISCARD constexpr line<3, ScalarT> intersection(hyperplane<3, ScalarT> cons
     return {p, dir};
 }
 
+template <int D, class ScalarT>
+TG_NODISCARD constexpr optional<aabb<D, ScalarT>> intersection(aabb<D, ScalarT> const& a, aabb<D, ScalarT> const& b)
+{
+    for (auto i = 0; i < D; ++i)
+    {
+        if (a.max[i] < b.min[i])
+            return {};
+
+        if (b.max[i] < a.min[i])
+            return {};
+    }
+
+    // overlap!
+    aabb<D, ScalarT> res;
+
+    for (auto i = 0; i < D; ++i)
+    {
+        res.min[i] = max(a.min[i], b.min[i]);
+        res.max[i] = min(a.max[i], b.max[i]);
+    }
+
+    return {res};
+}
+
 template <class ScalarT>
 TG_NODISCARD constexpr ScalarT intersection_parameter(line<3, ScalarT> const& l, hyperplane<3, ScalarT> const& p)
 {

src/typed-geometry/functions/norm.hh

@@ -33,14 +33,21 @@ TG_NODISCARD constexpr T norm(vec<4, T> const& v, T p)
 }
 
 template <int C, int R, class T>
-TG_NODISCARD constexpr T frobenius_norm(mat<C, R, T> const& v)
+TG_NODISCARD constexpr T frobenius_norm_sqr(mat<C, R, T> const& v)
 {
-    static_assert(is_floating_point<T>, "frobenius_norm only works on f32/f64");
+    static_assert(is_floating_point<T>, "frobenius_norm_sqr only works on f32/f64");
     auto s = T(0);
     for (auto x = 0; x < C; ++x)
         for (auto y = 0; y < R; ++y)
             s += tg::pow2(v[x][y]);
-    return tg::sqrt(s);
+    return s;
+}
+
+template <int C, int R, class T>
+TG_NODISCARD constexpr T frobenius_norm(mat<C, R, T> const& v)
+{
+    static_assert(is_floating_point<T>, "frobenius_norm only works on f32/f64");
+    return tg::sqrt(frobenius_norm_sqr(v));
 }
 
 }

src/typed-geometry/tests/triangle_tests.hh

+#pragma once
+
+#include <typed-geometry/detail/utility.hh>
+#include <typed-geometry/functions/distance.hh>
+#include <typed-geometry/types/objects/triangle.hh>
+
+namespace tg
+{
+enum class triangle_obtuseness
+{
+    obtuse, // largest angle > 90°
+    right,  // largest angle = 90°
+    acute   // largest angles < 90°
+};
+
+template <int D, class ScalarT>
+TG_NODISCARD constexpr triangle_obtuseness classify_obtuseness(triangle<D, ScalarT> const& t)
+{
+    auto e0sqr = tg::distance_sqr(t.pos0, t.pos1);
+    auto e1sqr = tg::distance_sqr(t.pos1, t.pos2);
+    auto e2sqr = tg::distance_sqr(t.pos2, t.pos0);
+
+    // make e0sqr longest (squared) edge
+    if (e1sqr > e0sqr)
+        tg::detail::swap(e0sqr, e1sqr);
+    if (e2sqr > e0sqr)
+        tg::detail::swap(e0sqr, e2sqr);
+
+    auto test = e0sqr - e1sqr - e2sqr;
+    if (test > ScalarT(0))
+        return triangle_obtuseness::obtuse;
+    if (test < ScalarT(0))
+        return triangle_obtuseness::acute;
+    return triangle_obtuseness::right;
+}
+
+template <int D, class ScalarT>
+TG_NODISCARD constexpr triangle_obtuseness is_obtuse_triangle(triangle<D, ScalarT> const& t)
+{
+    return classify_obtuseness(t) == triangle_obtuseness::obtuse;
+}
+
+template <int D, class ScalarT>
+TG_NODISCARD constexpr triangle_obtuseness is_right_triangle(triangle<D, ScalarT> const& t)
+{
+    return classify_obtuseness(t) == triangle_obtuseness::right;
+}
+
+template <int D, class ScalarT>
+TG_NODISCARD constexpr triangle_obtuseness is_nonobtuse_triangle(triangle<D, ScalarT> const& t)
+{
+    return classify_obtuseness(t) != triangle_obtuseness::obtuse;
+}
+
+template <int D, class ScalarT>
+TG_NODISCARD constexpr triangle_obtuseness is_acute_triangle(triangle<D, ScalarT> const& t)
+{
+    return classify_obtuseness(t) == triangle_obtuseness::acute;
+}
+}