Updated from develop

src/typed-geometry/detail/scalar_traits.hh

@@ -113,4 +113,20 @@ struct is_scalar_t<angle_t<T>>
 // abstract scalars are scalars that don't hold a concrete values (e.g. traced types)
 TG_IMPL_DEFINE_TRAIT(is_abstract_scalar, bool, false);
 
+template <class T>
+struct make_unsigned_t
+{
+    using type = std::make_unsigned_t<T>;
+};
+template <class T>
+using make_unsigned = typename make_unsigned_t<T>::type;
+
+template <class T>
+struct make_signed_t
+{
+    using type = std::make_signed_t<T>;
+};
+template <class T>
+using make_signed = typename make_signed_t<T>::type;
+
 } // namespace tg

src/typed-geometry/detail/utility.hh

@@ -253,6 +253,27 @@ struct is_range_t<Container,
                       >> : std::true_type
 {
 };
+
+template <class Container, class = void>
+struct is_any_range_t : false_type
+{
+};
+template <class ElementT, size_t N>
+struct is_any_range_t<ElementT[N]> : true_type
+{
+};
+template <class ElementT, size_t N>
+struct is_any_range_t<ElementT (&)[N]> : true_type
+{
+};
+template <class Container>
+struct is_any_range_t<Container,
+                      std::void_t<                                     //
+                          decltype(std::declval<Container>().begin()), //
+                          decltype(std::declval<Container>().end())    //
+                          >> : std::true_type
+{
+};
 }
 
 template <class Container, class ElementT>
@@ -261,6 +282,9 @@ static constexpr bool is_container = sizeof(detail::container_test<Container, El
 template <class Container, class ElementT>
 static constexpr bool is_range = detail::is_range_t<Container, ElementT>::value;
 
+template <class Container>
+static constexpr bool is_any_range = detail::is_any_range_t<Container>::value;
+
 template <class C>
 constexpr auto begin(C& c) -> decltype(c.begin())
 {

src/typed-geometry/feature/matrix.hh

@@ -3,7 +3,9 @@
 #include <typed-geometry/feature/vector.hh>
 
 #include <typed-geometry/functions/matrix/adjugate.hh>
+#include <typed-geometry/functions/matrix/compose.hh>
 #include <typed-geometry/functions/matrix/covariance.hh>
+#include <typed-geometry/functions/matrix/decompose.hh>
 #include <typed-geometry/functions/matrix/determinant.hh>
 #include <typed-geometry/functions/matrix/diag.hh>
 #include <typed-geometry/functions/matrix/eigenvalues.hh>

src/typed-geometry/feature/objects.hh

@@ -23,5 +23,6 @@
 #include <typed-geometry/functions/objects/rasterize.hh>
 #include <typed-geometry/functions/objects/size.hh>
 #include <typed-geometry/functions/objects/tangent.hh>
+#include <typed-geometry/functions/objects/triangle.hh>
 #include <typed-geometry/functions/objects/vertices.hh>
 #include <typed-geometry/functions/objects/volume.hh>

src/typed-geometry/functions/fixed_int/conversions.hh

 #pragma once
 
+#include <typed-geometry/detail/scalar_traits.hh>
 #include <typed-geometry/types/scalars/fixed_int.hh>
 #include <typed-geometry/types/scalars/fixed_uint.hh>
 
@@ -30,4 +31,45 @@ constexpr fixed_uint<words>::fixed_uint(fixed_int<rhs_words> const& rhs)
     if constexpr (rhs_words > 3 && words > 3)
         d[3] = rhs.d[3];
 }
+
+template <int w>
+struct make_unsigned_t<fixed_int<w>>
+{
+    using type = fixed_uint<w>;
+};
+template <int w>
+struct make_unsigned_t<fixed_int<w> const>
+{
+    using type = fixed_uint<w> const;
+};
+template <int w>
+struct make_unsigned_t<fixed_uint<w>>
+{
+    using type = fixed_uint<w>;
+};
+template <int w>
+struct make_unsigned_t<fixed_uint<w> const>
+{
+    using type = fixed_uint<w> const;
+};
+template <int w>
+struct make_signed_t<fixed_uint<w>>
+{
+    using type = fixed_int<w>;
+};
+template <int w>
+struct make_signed_t<fixed_uint<w> const>
+{
+    using type = fixed_int<w> const;
+};
+template <int w>
+struct make_signed_t<fixed_int<w>>
+{
+    using type = fixed_int<w>;
+};
+template <int w>
+struct make_signed_t<fixed_int<w> const>
+{
+    using type = fixed_int<w> const;
+};
 }

src/typed-geometry/functions/matrix/compose.hh

+#pragma once
+
+#include <typed-geometry/detail/operators.hh>
+#include <typed-geometry/functions/basic/scalar_math.hh>
+#include <typed-geometry/types/mat.hh>
+#include <typed-geometry/types/quat.hh>
+#include <typed-geometry/types/size.hh>
+#include <typed-geometry/types/vec.hh>
+
+namespace tg
+{
+template <class T>
+constexpr mat<4, 4, T> compose_transformation(vec<3, T> const& translation, mat<3, 3, T> const& rotation, float scale)
+{
+    mat<4, 4, T> M;
+    M[0] = tg::vec4(rotation[0] * scale, T(0));
+    M[1] = tg::vec4(rotation[1] * scale, T(0));
+    M[2] = tg::vec4(rotation[2] * scale, T(0));
+    M[3] = tg::vec4(translation, T(1));
+    return M;
+}
+template <class T>
+constexpr mat<4, 4, T> compose_transformation(vec<3, T> const& translation, mat<3, 3, T> const& rotation, size<3, T> const& scale)
+{
+    mat<4, 4, T> M;
+    M[0] = tg::vec4(rotation[0] * scale[0], T(0));
+    M[1] = tg::vec4(rotation[1] * scale[1], T(0));
+    M[2] = tg::vec4(rotation[2] * scale[2], T(0));
+    M[3] = tg::vec4(translation, T(1));
+    return M;
+}
+}

src/typed-geometry/functions/matrix/decompose.hh

+#pragma once
+
+#include <typed-geometry/detail/operators.hh>
+#include <typed-geometry/functions/basic/scalar_math.hh>
+#include <typed-geometry/types/mat.hh>
+#include <typed-geometry/types/quat.hh>
+#include <typed-geometry/types/size.hh>
+#include <typed-geometry/types/vec.hh>
+
+namespace tg
+{
+/// decomposes a transformation matrix into a translation t, rotation matrix R, and uniform scale s
+///
+/// Satisfies:
+///   M = tg::translation(t) * tg::rotation(R) * tg::scaling<3>(s)
+///   or:
+///   M * p = R * s * p + t
+///
+/// Usage:
+///
+///   tg::mat4 M = ...;
+///   tg::vec3 t;
+///   tg::mat3 R;
+///   float s;
+///   decompose_transformation(M, t, R, s);
+///   .. use t,R,s
+///
+/// NOTE: assumes that M is actually properly decomposable
+///       (e.g. rotation may not be orthonormal if M has non-uniform scale)
+template <class T>
+constexpr void decompose_transformation(mat<4, 4, T> const& M, vec<3, T>& translation, mat<3, 3, T>& rotation, float& scale)
+{
+    // last column is translation
+    translation[0] = M[3][0];
+    translation[1] = M[3][1];
+    translation[2] = M[3][2];
+
+    // uniform scale is length of first col
+    auto s = length(vec<3, T>(M[0]));
+    scale = s;
+
+    // rot is 3x3 part
+    auto inv_s = T(1) / s;
+    rotation[0] = vec<3, T>(M[0]) * inv_s;
+    rotation[1] = vec<3, T>(M[1]) * inv_s;
+    rotation[2] = vec<3, T>(M[2]) * inv_s;
+}
+template <class T>
+constexpr void decompose_transformation(mat<4, 4, T> const& M, vec<3, T>& translation, mat<3, 3, T>& rotation, size<3, T>& scale)
+{
+    // last column is translation
+    translation[0] = M[3][0];
+    translation[1] = M[3][1];
+    translation[2] = M[3][2];
+
+    // scale is length of cols
+    auto s0 = length(vec<3, T>(M[0]));
+    auto s1 = length(vec<3, T>(M[1]));
+    auto s2 = length(vec<3, T>(M[2]));
+    scale = {s0, s1, s2};
+
+    // rot is 3x3 part
+    rotation[0] = vec<3, T>(M[0]) / s0;
+    rotation[1] = vec<3, T>(M[1]) / s1;
+    rotation[2] = vec<3, T>(M[2]) / s2;
+}
+}

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

@@ -102,4 +102,16 @@ template <class ObjectT>
 {
     return edges_of(o);
 }
+
+namespace detail
+{
+template <class T>
+auto test_has_edges(int) -> decltype(void(edges_of(std::declval<T>())), std::true_type{});
+template <class T>
+std::false_type test_has_edges(char);
+}
+
+/// true if edges(obj) exists
+template <class T>
+constexpr bool has_edges_of = decltype(detail::test_has_edges<T>(0))::value;
 }

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

@@ -83,4 +83,16 @@ template <class BaseT, class TraitsT>
     auto mantle = faces_of(boundary_no_caps_of(py));
     return pyramid_faces<BaseT, mantle.size()>(py.base, mantle);
 }
+
+namespace detail
+{
+template <class T>
+auto test_has_faces(int) -> decltype(void(faces_of(std::declval<T>())), std::true_type{});
+template <class T>
+std::false_type test_has_faces(char);
+}
+
+/// true if faces_of(obj) exists
+template <class T>
+constexpr bool has_faces_of = decltype(detail::test_has_faces<T>(0))::value;
 }

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

 #pragma once
 
+#include <typed-geometry/functions/matrix/inverse.hh>
 #include <typed-geometry/functions/objects/distance.hh>
+#include <typed-geometry/functions/objects/normal.hh>
 #include <typed-geometry/types/objects/sphere.hh>
 #include <typed-geometry/types/objects/triangle.hh>
 
 namespace tg
 {
-template <class ScalarT>
-[[nodiscard]] sphere<2, ScalarT> circumcircle_of(triangle<2, ScalarT> const& t)
+template <int D, class ScalarT>
+[[nodiscard]] pos<D, ScalarT> circumcenter_of(triangle<D, ScalarT> const& t)
 {
-    auto const& a = t.pos0;
-    auto const& b = t.pos1;
-    auto const& c = t.pos2;
-
-    auto const d = ScalarT(2) * (a.x * (b.y - c.y) + b.x * (c.y - a.y) + c.x * (a.y - b.y));
-    auto const a_len_sqr = distance_sqr_to_origin(a);
-    auto const b_len_sqr = distance_sqr_to_origin(b);
-    auto const c_len_sqr = distance_sqr_to_origin(c);
-    auto const x = (a_len_sqr * (b.y - c.y) + b_len_sqr * (c.y - a.y) + c_len_sqr * (a.y - b.y)) / d;
-    auto const y = (a_len_sqr * (c.x - b.x) + b_len_sqr * (a.x - c.x) + c_len_sqr * (b.x - a.x)) / d;
+    auto const a2 = distance_sqr(t.pos0, t.pos1);
+    auto const b2 = distance_sqr(t.pos1, t.pos2);
+    auto const c2 = distance_sqr(t.pos2, t.pos0);
 
-    auto const center = pos<2, ScalarT>(x, y);
-    auto const r = distance(center, a);
+    auto A = a2 * (b2 + c2 - a2);
+    auto B = b2 * (c2 + a2 - b2);
+    auto C = c2 * (a2 + b2 - c2);
+    auto const f = ScalarT(1) / (A + B + C);
+    A *= f;
+    B *= f;
+    C *= f;
+    return t.pos0 * B + t.pos1 * C + t.pos2 * A;
+}
 
-    return {center, r};
+template <int D, class ScalarT>
+[[nodiscard]] ScalarT circumradius_of(triangle<D, ScalarT> const& t)
+{
+    auto const a2 = distance_sqr(t.pos0, t.pos1);
+    auto const b2 = distance_sqr(t.pos1, t.pos2);
+    auto const c2 = distance_sqr(t.pos2, t.pos0);
+    return sqrt(a2 * b2 * c2) / (4 * area_of(t));
 }
 
+template <int D, class ScalarT>
+[[nodiscard]] sphere<2, ScalarT, D> circumcircle_of(triangle<D, ScalarT> const& t)
+{
+    auto const center = circumcenter_of(t);
+    auto const r = distance(center, t.pos0);
+
+    if constexpr (D == 2)
+        return {center, r};
+    else
+        return {center, r, normal_of(t)};
+}
 
 template <class ObjectT>
 [[deprecated("use circumcircle_of")]] [[nodiscard]] constexpr auto circumcircle(ObjectT const& o) -> decltype(circumcircle_of(o))
 {
     return circumcircle_of(o);
 }
+
+template <int D, class ScalarT>
+[[nodiscard]] pos<D, ScalarT> incenter_of(triangle<D, ScalarT> const& t)
+{
+    auto a = distance(t.pos0, t.pos1);
+    auto b = distance(t.pos1, t.pos2);
+    auto c = distance(t.pos2, t.pos0);
+    auto const f = ScalarT(1) / (a + b + c);
+    a *= f;
+    b *= f;
+    c *= f;
+    return t.pos0 * b + t.pos1 * c + t.pos2 * a;
+}
+
+template <int D, class ScalarT>
+[[nodiscard]] ScalarT inradius_of(triangle<D, ScalarT> const& t)
+{
+    auto const a = distance(t.pos0, t.pos1);
+    auto const b = distance(t.pos1, t.pos2);
+    auto const c = distance(t.pos2, t.pos0);
+    auto const s = (a + b + c) / 2;
+    return sqrt(max(ScalarT(0), (s - a) * (s - b) * (s - c) / s));
+}
+
+template <int D, class ScalarT>
+[[nodiscard]] sphere<2, ScalarT, D> incircle_of(triangle<D, ScalarT> const& t)
+{
+    auto const center = incenter_of(t);
+    auto const r = inradius_of(t);
+
+    if constexpr (D == 2)
+        return {center, r};
+    else
+        return {center, r, normal_of(t)};
+}
+
+/// given a 2D triangle, returns a 3x3 matrix M such that
+/// M * (a, b, c) = (p, 1)
+/// where a,b,c are barycentric coordinates for p
+template <class ScalarT>
+[[nodiscard]] mat<3, 3, ScalarT> from_barycoord_matrix_of(triangle<2, ScalarT> const& t)
+{
+    return mat<3, 3, ScalarT>::from_cols(vec<3, ScalarT>(t.pos0, 1.f), //
+                                         vec<3, ScalarT>(t.pos1, 1.f), //
+                                         vec<3, ScalarT>(t.pos2, 1.f));
+}
+
+/// given a 3D triangle, returns a 4x4 matrix M such that
+/// M * (a, b, c, d) = (p + normal * d, 1)
+/// where a,b,c are barycentric coordinates for p
+/// and d is the signed distance in normal direction
+template <class ScalarT>
+[[nodiscard]] mat<4, 4, ScalarT> from_barycoord_matrix_of(triangle<3, ScalarT> const& t)
+{
+    auto const n = normal_of(t);
+    return mat<4, 4, ScalarT>::from_cols(vec<4, ScalarT>(t.pos0, 1.f), //
+                                         vec<4, ScalarT>(t.pos1, 1.f), //
+                                         vec<4, ScalarT>(t.pos2, 1.f), //
+                                         vec<4, ScalarT>(n, 0.f));
+}
+
+/// given a 2D triangle, returns a 3x3 matrix M such that
+/// M * (p, 1) = (a, b, c)
+/// where a,b,c are barycentric coordinates for p
+template <class ScalarT>
+[[nodiscard]] mat<3, 3, ScalarT> to_barycoord_matrix_of(triangle<2, ScalarT> const& t)
+{
+    return inverse(from_barycoord_matrix_of(t));
+}
+
+/// given a 3D triangle, returns a 4x4 matrix M such that
+/// M * (p, 1) = (a, b, c, d)
+/// where a,b,c are barycentric coordinates for p
+/// and d is the signed distance in normal direction
+template <class ScalarT>
+[[nodiscard]] mat<4, 4, ScalarT> to_barycoord_matrix_of(triangle<3, ScalarT> const& t)
+{
+    return inverse(from_barycoord_matrix_of(t));
+}
+
+/// returns a transformation matrix M such that
+/// M * p = p'
+/// where p is a point in "from" and p' in "to" with the same barycentric coordinates
+/// (i.e. it's the matrix transforming the first triangle into the second)
+template <int D, class ScalarT>
+[[nodiscard]] mat<D + 1, D + 1, ScalarT> transformation_from_to(triangle<D, ScalarT> const& from, triangle<D, ScalarT> const& to)
+{
+    return from_barycoord_matrix_of(to) * to_barycoord_matrix_of(from);
+}
+
+/// returns the minimum altitude/height of a triangle
+template <int D, class ScalarT>
+[[nodiscard]] ScalarT min_height_of(triangle<D, ScalarT> const& t)
+{
+    auto const area = area_of(t);
+    return area / tg::max(distance(t.pos0, t.pos1), distance(t.pos1, t.pos2), distance(t.pos2, t.pos0));
+}
+
+/// returns the maximum altitude/height of a triangle
+template <int D, class ScalarT>
+[[nodiscard]] ScalarT max_height_of(triangle<D, ScalarT> const& t)
+{
+    auto const area = area_of(t);
+    return area / tg::min(distance(t.pos0, t.pos1), distance(t.pos1, t.pos2), distance(t.pos2, t.pos0));
+}
 }

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

@@ -100,4 +100,16 @@ template <class ObjectT>
 {
     return vertices_of(o);
 }
+
+namespace detail
+{
+template <class T>
+auto test_has_vertices(int) -> decltype(void(vertices_of(std::declval<T>())), std::true_type{});
+template <class T>
+std::false_type test_has_vertices(char);
+}
+
+/// true if vertices_of(obj) exists
+template <class T>
+constexpr bool has_vertices_of = decltype(detail::test_has_vertices<T>(0))::value;
 }

src/typed-geometry/functions/random/uniform.hh

@@ -748,6 +748,20 @@ template <class Rng, class... Args>
     return uniform(rng, args...) - decltype(uniform(rng, args...))::zero;
 }
 
+/// returns uniformly sampled barycentric coordinates
+/// i.e. uniform(rng, triangle) has the same distribution as
+///      triangle[uniform_barycoords(rng)]
+template <class ScalarT = float, class Rng>
+[[nodiscard]] constexpr comp<3, ScalarT> uniform_barycoords(Rng& rng)
+{
+    auto const a = detail::uniform01<ScalarT>(rng);
+    auto const b = detail::uniform01<ScalarT>(rng);
+    if (a + b <= ScalarT(1))
+        return {a, b, 1 - a - b};
+    else
+        return {1 - a, 1 - b, a + b - 1};
+}
+
 
 namespace detail
 {

src/typed-geometry/functions/std/io.hh

@@ -73,37 +73,31 @@ std::basic_ostream<CharT, StreamTraits>& operator<<(std::basic_ostream<CharT, St
 template <int w, class CharT, class StreamTraits>
 std::basic_ostream<CharT, StreamTraits>& operator<<(std::basic_ostream<CharT, StreamTraits>& out, fixed_uint<w> const& val)
 {
-    constexpr const char* hex_map = "0123456789ABCDEF";
-
-    auto ss = detail::temp_sstream(out);
-    ss << "u" << 64 * w << "(";
-    ss << "0x";
-    for (auto wi = w - 1; wi >= 0; --wi)
-        for (auto i = 15; i >= 0; --i)
-        {
-            const u64 idx = (val.d[wi] >> (i * 4)) & 0xF;
-            ss << hex_map[idx];
-        }
-    ss << ")";
-    return out << ss.str();
+    CharT buf[w * 20 + 1];
+    auto const end = buf + sizeof(buf);
+    auto begin = end;
+    *(--begin) = '\0';
+    auto u = val;
+    while (u > 0)
+    {
+        auto const n = u64(u % 10);
+        u /= 10;
+        *(--begin) = CharT('0' + n);
+    }
+    if (begin + 1 == end)
+        *(--begin) = CharT('0');
+    return out << begin;
 }
 
 template <int w, class CharT, class StreamTraits>
 std::basic_ostream<CharT, StreamTraits>& operator<<(std::basic_ostream<CharT, StreamTraits>& out, fixed_int<w> const& val)
 {
-    constexpr const char* hex_map = "0123456789ABCDEF";
-
-    auto ss = detail::temp_sstream(out);
-    ss << "i" << 64 * w << "(";
-    ss << "0x";
-    for (auto wi = w - 1; wi >= 0; --wi)
-        for (auto i = 15; i >= 0; --i)
-        {
-            const u64 idx = (val.d[wi] >> (i * 4)) & 0xF;
-            ss << hex_map[idx];
-        }
-    ss << ")";
-    return out << ss.str();
+    if (val < 0)
+    {
+        out << '-';
+        return out << fixed_uint<w>(-val);
+    }
+    return out << fixed_uint<w>(val);
 }
 
 //