From e06765346ebd8dc2f7130c7e8bf6f8f29f07db1f Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Tue, 15 Oct 2024 11:22:03 -0700 Subject: [PATCH] Replace GLM with linalg (#984) * added linalg * removed GLM * fix matrix multiplication * more matrix constructors * fixed more functions * fixed matrix notation * more fixes * fix ostream * further fixes * added more rotation functions * compiles * some fixes * fix2 * fix build * another warning... * fix install * for wasm? * no warning for unknown options * proper fix? * remove legacy warning flags * format * misc constexpr --------- Co-authored-by: pca006132 --- CMakeLists.txt | 28 +- README.md | 1 - bindings/c/box.cpp | 2 +- bindings/c/cross.cpp | 2 +- bindings/c/manifoldc.cpp | 2 +- bindings/c/rect.cpp | 2 +- bindings/python/manifold3d.cpp | 64 +- bindings/wasm/helpers.cpp | 4 +- extras/minimize_testcase.cpp | 10 +- extras/test_hull_performance.cpp | 5 +- flake.nix | 3 - include/manifold/common.h | 184 +-- include/manifold/cross_section.h | 4 +- include/manifold/iters.h | 75 +- include/manifold/linalg.h | 2095 +++++++++++++++++++++++++++ include/manifold/manifold.h | 20 +- manifold.pc.in | 2 +- manifoldConfig.cmake.in | 2 - samples/src/bracelet.cpp | 14 +- samples/src/condensed_matter.cpp | 6 +- samples/src/gyroid_module.cpp | 8 +- samples/src/knot.cpp | 8 +- samples/src/scallop.cpp | 8 +- src/CMakeLists.txt | 1 - src/boolean3.cpp | 14 +- src/boolean_result.cpp | 8 +- src/collider.h | 6 +- src/constructors.cpp | 41 +- src/cross_section/cross_section.cpp | 38 +- src/csg_tree.cpp | 55 +- src/csg_tree.h | 18 +- src/edge_op.cpp | 18 +- src/face_op.cpp | 18 +- src/impl.cpp | 54 +- src/impl.h | 17 +- src/manifold.cpp | 34 +- src/meshIO/meshIO.cpp | 10 +- src/mesh_fixes.h | 2 +- src/polygon.cpp | 38 +- src/properties.cpp | 36 +- src/quickhull.cpp | 12 +- src/quickhull.h | 6 +- src/sdf.cpp | 30 +- src/shared.h | 48 +- src/smoothing.cpp | 161 +- src/subdivision.cpp | 48 +- src/svd.h | 24 +- src/tri_dist.h | 71 +- src/utils.h | 70 +- test/CMakeLists.txt | 13 +- test/boolean_complex_test.cpp | 30 +- test/cross_section_test.cpp | 24 +- test/hull_test.cpp | 4 +- test/manifold_test.cpp | 43 +- test/manifoldc_test.cpp | 9 +- test/polygon_test.cpp | 2 +- test/samples_test.cpp | 2 +- test/sdf_test.cpp | 38 +- test/smooth_test.cpp | 45 +- test/test.h | 7 - test/test_main.cpp | 44 +- 61 files changed, 2882 insertions(+), 806 deletions(-) create mode 100644 include/manifold/linalg.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 67aa27215..2607e79f2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -152,17 +152,16 @@ if(MSVC) list(APPEND MANIFOLD_FLAGS /DNOMINMAX /bigobj) else() if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU") + list(APPEND WARNING_FLAGS -Wall -Wno-unknown-warning-option -Wno-unused) + else() list( APPEND WARNING_FLAGS -Wall + -Wno-unknown-warning-option -Wno-unused - -Wno-array-bounds - -Wno-stringop-overflow - -Wno-alloc-size-larger-than + -Wno-shorten-64-to-32 ) - else() - list(APPEND WARNING_FLAGS -Wall -Wno-unused -Wno-shorten-64-to-32) endif() if(CMAKE_SYSTEM_NAME STREQUAL "Windows") if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU") @@ -218,24 +217,6 @@ function(logmissingdep PKG) endif() endfunction() -# GLM is required in all configurations. -find_package(glm QUIET) -if(NOT glm_FOUND) - logmissingdep("glm") - message(STATUS "glm not found, downloading from source") - set(GLM_BUILD_INSTALL "ON" CACHE STRING "") - FetchContent_Declare( - glm - GIT_REPOSITORY https://github.com/g-truc/glm.git - GIT_TAG 1.0.1 - GIT_PROGRESS TRUE - ) - FetchContent_MakeAvailable(glm) - if(NOT EMSCRIPTEN) - install(TARGETS glm) - endif() -endif() - # If we're building parallel, we need the requisite libraries if(MANIFOLD_PAR) find_package(Threads REQUIRED) @@ -411,6 +392,7 @@ set( MANIFOLD_PUBLIC_HDRS include/manifold/common.h include/manifold/iters.h + include/manifold/linalg.h include/manifold/manifold.h include/manifold/optional_assert.h include/manifold/parallel.h diff --git a/README.md b/README.md index b0adab2c4..10749d400 100644 --- a/README.md +++ b/README.md @@ -33,7 +33,6 @@ If you prefer Python to JS/TS, make your own copy of the example notebook above. This is a modern C++ library that Github's CI verifies builds and runs on a variety of platforms. Additionally, we build bindings for JavaScript ([manifold-3d](https://www.npmjs.com/package/manifold-3d) on npm), Python ([manifold3d](https://pypi.org/project/manifold3d/)), and C to make this library more portable and easy to use. System Dependencies (note that we will automatically download the dependency if there is no such package on the system): -- [`GLM`](https://github.com/g-truc/glm/): A compact header-only vector library. - [`tbb`](https://github.com/oneapi-src/oneTBB/): Intel's thread building blocks library. (only when `MANIFOLD_PAR=ON` is enabled) - [`gtest`](https://github.com/google/googletest/): Google test library (only when test is enabled, i.e. `MANIFOLD_TEST=ON`) diff --git a/bindings/c/box.cpp b/bindings/c/box.cpp index 81c6d0126..e471edba6 100644 --- a/bindings/c/box.cpp +++ b/bindings/c/box.cpp @@ -69,7 +69,7 @@ ManifoldBox *manifold_box_transform(void *mem, ManifoldBox *b, double x1, double y1, double z1, double x2, double y2, double z2, double x3, double y3, double z3, double x4, double y4, double z4) { - auto mat = mat4x3(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4); + auto mat = mat3x4({x1, y1, z1}, {x2, y2, z2}, {x3, y3, z3}, {x4, y4, z4}); auto transformed = from_c(b)->Transform(mat); return to_c(new (mem) Box(transformed)); } diff --git a/bindings/c/cross.cpp b/bindings/c/cross.cpp index 3142c80ef..d05e07e46 100644 --- a/bindings/c/cross.cpp +++ b/bindings/c/cross.cpp @@ -176,7 +176,7 @@ ManifoldCrossSection *manifold_cross_section_transform(void *mem, double x1, double y1, double x2, double y2, double x3, double y3) { - auto mat = mat3x2(x1, y1, x2, y2, x3, y3); + auto mat = mat2x3({x1, y1}, {x2, y2}, {x3, y3}); auto transformed = from_c(cs)->Transform(mat); return to_c(new (mem) CrossSection(transformed)); } diff --git a/bindings/c/manifoldc.cpp b/bindings/c/manifoldc.cpp index 1569bc061..20930cc20 100644 --- a/bindings/c/manifoldc.cpp +++ b/bindings/c/manifoldc.cpp @@ -237,7 +237,7 @@ ManifoldManifold *manifold_transform(void *mem, ManifoldManifold *m, double x1, double y1, double z1, double x2, double y2, double z2, double x3, double y3, double z3, double x4, double y4, double z4) { - auto mat = mat4x3(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4); + auto mat = mat3x4({x1, y1, z1}, {x2, y2, z2}, {x3, y3, z3}, {x4, y4, z4}); auto transformed = from_c(m)->Transform(mat); return to_c(new (mem) Manifold(transformed)); } diff --git a/bindings/c/rect.cpp b/bindings/c/rect.cpp index ff8d66fb9..86fe4ed8e 100644 --- a/bindings/c/rect.cpp +++ b/bindings/c/rect.cpp @@ -73,7 +73,7 @@ ManifoldRect *manifold_rect_union(void *mem, ManifoldRect *a, ManifoldRect *b) { ManifoldRect *manifold_rect_transform(void *mem, ManifoldRect *r, double x1, double y1, double x2, double y2, double x3, double y3) { - auto mat = mat3x2(x1, y1, x2, y2, x3, y3); + auto mat = mat2x3({x1, y1}, {x2, y2}, {x3, y3}); auto transformed = from_c(r)->Transform(mat); return to_c(new (mem) Rect(transformed)); } diff --git a/bindings/python/manifold3d.cpp b/bindings/python/manifold3d.cpp index 3eae563e5..9aaab4815 100644 --- a/bindings/python/manifold3d.cpp +++ b/bindings/python/manifold3d.cpp @@ -32,36 +32,36 @@ namespace nb = nanobind; using namespace manifold; template -struct glm_name {}; +struct la_name {}; template <> -struct glm_name { +struct la_name { static constexpr char const name[] = "Doublex3"; static constexpr char const multi_name[] = "DoubleNx3"; }; template <> -struct glm_name { +struct la_name { static constexpr char const name[] = "Doublex2"; static constexpr char const multi_name[] = "DoubleNx2"; }; template <> -struct glm_name { +struct la_name { static constexpr char const name[] = "Intx3"; static constexpr char const multi_name[] = "IntNx3"; }; template <> -struct glm_name { +struct la_name { static constexpr char const name[] = "Double3x4"; }; template <> -struct glm_name { +struct la_name { static constexpr char const name[] = "Double2x3"; }; -// handle glm::vecN -template -struct nb::detail::type_caster> { - using glm_type = glm::vec; - NB_TYPE_CASTER(glm_type, const_name(glm_name::name)); +// handle la::vecN +template +struct nb::detail::type_caster> { + using la_type = la::vec; + NB_TYPE_CASTER(la_type, const_name(la_name::name)); bool from_python(handle src, uint8_t flags, cleanup_list *cleanup) noexcept { int size = PyObject_Size(src.ptr()); // negative on failure @@ -73,7 +73,7 @@ struct nb::detail::type_caster> { } return true; } - static handle from_cpp(glm_type vec, rv_policy policy, + static handle from_cpp(la_type vec, rv_policy policy, cleanup_list *cleanup) noexcept { nb::list out; for (int i = 0; i < N; i++) out.append(vec[i]); @@ -81,12 +81,12 @@ struct nb::detail::type_caster> { } }; -// handle glm::matMxN -template -struct nb::detail::type_caster> { - using glm_type = glm::mat; +// handle la::matMxN +template +struct nb::detail::type_caster> { + using la_type = la::mat; using numpy_type = nb::ndarray>; - NB_TYPE_CASTER(glm_type, const_name(glm_name::name)); + NB_TYPE_CASTER(la_type, const_name(la_name::name)); bool from_python(handle src, uint8_t flags, cleanup_list *cleanup) noexcept { int rows = PyObject_Size(src.ptr()); // negative on failure @@ -103,13 +103,13 @@ struct nb::detail::type_caster> { } return true; } - static handle from_cpp(glm_type mat, rv_policy policy, + static handle from_cpp(la_type mat, rv_policy policy, cleanup_list *cleanup) noexcept { T *buffer = new T[R * C]; nb::capsule mem_mgr(buffer, [](void *p) noexcept { delete[] (T *)p; }); for (int i = 0; i < R; i++) { for (int j = 0; j < C; j++) { - // py is (Rows, Cols), glm is (Cols, Rows) + // py is (Rows, Cols), la is (Cols, Rows) buffer[i * C + j] = mat[j][i]; } } @@ -119,13 +119,13 @@ struct nb::detail::type_caster> { } }; -// handle std::vector -template -struct nb::detail::type_caster>> { - using glm_type = glm::vec; +// handle std::vector +template +struct nb::detail::type_caster>> { + using la_type = la::vec; using numpy_type = nb::ndarray>; - NB_TYPE_CASTER(std::vector, - const_name(glm_name::multi_name)); + NB_TYPE_CASTER(std::vector, + const_name(la_name::multi_name)); bool from_python(handle src, uint8_t flags, cleanup_list *cleanup) noexcept { make_caster arr_cast; @@ -142,7 +142,7 @@ struct nb::detail::type_caster>> { if (num_vec == static_cast(-1)) return false; value.resize(num_vec); for (size_t i = 0; i < num_vec; i++) { - make_caster vec_cast; + make_caster vec_cast; if (!vec_cast.from_python(src[i], flags, cleanup)) return false; value[i] = vec_cast.value; } @@ -165,13 +165,13 @@ struct nb::detail::type_caster>> { } }; -// handle VecView -template -struct nb::detail::type_caster>> { - using glm_type = glm::vec; +// handle VecView +template +struct nb::detail::type_caster>> { + using la_type = la::vec; using numpy_type = nb::ndarray>; - NB_TYPE_CASTER(manifold::VecView, - const_name(glm_name::multi_name)); + NB_TYPE_CASTER(manifold::VecView, + const_name(la_name::multi_name)); bool from_python(handle src, uint8_t flags, cleanup_list *cleanup) noexcept { make_caster arr_cast; diff --git a/bindings/wasm/helpers.cpp b/bindings/wasm/helpers.cpp index 8b2c82bac..3e442c6fe 100644 --- a/bindings/wasm/helpers.cpp +++ b/bindings/wasm/helpers.cpp @@ -140,7 +140,7 @@ CrossSection IntersectionN(const std::vector& cross_sections) { CrossSection Transform(CrossSection& cross_section, const val& mat) { std::vector array = convertJSArrayToNumberVector(mat); - mat3x2 matrix; + mat2x3 matrix; for (const int col : {0, 1, 2}) for (const int row : {0, 1}) matrix[col][row] = array[col * 3 + row]; return cross_section.Transform(matrix); @@ -188,7 +188,7 @@ Manifold IntersectionN(const std::vector& manifolds) { Manifold Transform(Manifold& manifold, const val& mat) { std::vector array = convertJSArrayToNumberVector(mat); - mat4x3 matrix; + mat3x4 matrix; for (const int col : {0, 1, 2, 3}) for (const int row : {0, 1, 2}) matrix[col][row] = array[col * 4 + row]; return manifold.Transform(matrix); diff --git a/extras/minimize_testcase.cpp b/extras/minimize_testcase.cpp index 9242b61a5..85561d3b8 100644 --- a/extras/minimize_testcase.cpp +++ b/extras/minimize_testcase.cpp @@ -23,8 +23,8 @@ inline bool intersect(vec2 p0, vec2 p1, vec2 q0, vec2 q1, double precision) { vec2 r = p1 - p0; vec2 s = q1 - q0; // allow some error in the intersection point - double epsilon_r = 2.0 * precision / glm::length(r); - double epsilon_s = 2.0 * precision / glm::length(s); + double epsilon_r = 2.0 * precision / la::length(r); + double epsilon_s = 2.0 * precision / la::length(s); double rxs = cross(r, s); // in case they are nearly collinear, ignore them... // this is to avoid treating degenerate triangles as intersecting @@ -45,7 +45,7 @@ inline bool intersect(vec2 p0, vec2 p1, vec2 q0, vec2 q1, double precision) { // perturbation along r/s is not enough. // // in that case, apply perturbation perpendicular to r - vec2 r_orth = glm::normalize(vec2(-r.y, r.x)) * precision; + vec2 r_orth = la::normalize(vec2(-r.y, r.x)) * precision; double u1 = cross(q0 - p0 + r_orth, r) / rxs; double u2 = cross(q0 - p0 - r_orth, r) / rxs; double t1 = cross(q0 - p0 + r_orth, s) / rxs; @@ -369,8 +369,8 @@ int main(int argc, char **argv) { double bound = 0; for (const SimplePolygon &poly : polys) { for (const vec2 &pt : poly) { - bound = glm::max(bound, glm::abs(pt.x)); - bound = glm::max(bound, glm::abs(pt.y)); + bound = la::max(bound, la::abs(pt.x)); + bound = la::max(bound, la::abs(pt.y)); } } precision = bound * kTolerance; diff --git a/extras/test_hull_performance.cpp b/extras/test_hull_performance.cpp index d4d2d522d..8df799576 100644 --- a/extras/test_hull_performance.cpp +++ b/extras/test_hull_performance.cpp @@ -31,7 +31,6 @@ #include "manifold/meshIO.h" #include "samples.h" using namespace std; -using namespace glm; using namespace manifold; // Epick = Exact predicates Inexact constructions. Seems fair to use to compare @@ -70,7 +69,7 @@ class HullImpl { vec3 v2 = vertPos[tri[2]]; // Compute the normal of the triangle - vec3 normal = glm::normalize(glm::cross(v1 - v0, v2 - v0)); + vec3 normal = la::normalize(la::cross(v1 - v0, v2 - v0)); // Check all other vertices for (int i = 0; i < (int)vertPos.size(); ++i) { @@ -81,7 +80,7 @@ class HullImpl { vec3 v = vertPos[i]; // Compute the signed distance from the plane - double distance = glm::dot(normal, v - v0); + double distance = la::dot(normal, v - v0); // If any vertex lies on the opposite side of the normal direction if (distance > 0) { diff --git a/flake.nix b/flake.nix index 02d27d328..2a1a8a5e9 100644 --- a/flake.nix +++ b/flake.nix @@ -40,7 +40,6 @@ pkg-config ]) ++ build-tools; buildInputs = with pkgs; [ - glm clipper2 assimp ]; @@ -86,7 +85,6 @@ mkdir build cd build emcmake cmake -DCMAKE_BUILD_TYPE=Release \ - -DFETCHCONTENT_SOURCE_DIR_GLM=${pkgs.glm.src} \ -DFETCHCONTENT_SOURCE_DIR_GOOGLETEST=${gtest-src} \ -DFETCHCONTENT_SOURCE_DIR_CLIPPER2=../clipper2 .. ''; @@ -114,7 +112,6 @@ ]; buildInputs = with pkgs; [ tbb - glm clipper2 ]; nativeBuildInputs = with pkgs; [ diff --git a/include/manifold/common.h b/include/manifold/common.h index ec3e553c3..7bcdc486a 100644 --- a/include/manifold/common.h +++ b/include/manifold/common.h @@ -13,59 +13,71 @@ // limitations under the License. #pragma once -#define GLM_ENABLE_EXPERIMENTAL // needed for glm/gtx/compatibility.hpp -#define GLM_FORCE_EXPLICIT_CTOR -#include -#include -#include -#include -#include #include #include #include +#include "manifold/linalg.h" + namespace manifold { /** @defgroup Math data structure definitions - * @brief Abstract away from glm. + * @brief Abstract away from linalg. * In the future the underlying data type can change. * @{ */ -using vec2 = glm::dvec2; -using vec3 = glm::dvec3; -using vec4 = glm::dvec4; -using mat2 = glm::dmat2; -using mat2x3 = glm::dmat2x3; -using mat2x4 = glm::dmat2x4; -using mat3x2 = glm::dmat3x2; -using mat3 = glm::dmat3; -using mat3x4 = glm::dmat3x4; -using mat4x3 = glm::dmat4x3; -using mat4 = glm::dmat4; -using ivec2 = glm::vec<2, int>; -using ivec3 = glm::vec<3, int>; -using ivec4 = glm::vec<4, int>; -using quat = glm::dquat; +namespace la = linalg; +using vec2 = la::vec; +using vec3 = la::vec; +using vec4 = la::vec; +using bvec4 = la::vec; +using mat2 = la::mat; +using mat3x2 = la::mat; +using mat4x2 = la::mat; +using mat2x3 = la::mat; +using mat3 = la::mat; +using mat4x3 = la::mat; +using mat3x4 = la::mat; +using mat4 = la::mat; +using ivec2 = la::vec; +using ivec3 = la::vec; +using ivec4 = la::vec; +using quat = la::vec; ///@} +constexpr double kPi = 3.14159265358979323846264338327950288; +constexpr double kTwoPi = 6.28318530717958647692528676655900576; +constexpr double kHalfPi = 1.57079632679489661923132169163975144; + +constexpr double radians(double a) { return a * kPi / 180; } +constexpr double degrees(double a) { return a * 180 / kPi; } + +constexpr double smoothstep(double edge0, double edge1, double a) { + const double x = la::clamp((a - edge0) / (edge1 - edge0), 0, 1); + return x * x * (3 - 2 * x); +} + +constexpr mat3x4 Identity3x4() { return mat3x4(mat3(la::identity), vec3(0.0)); } +constexpr mat2x3 Identity2x3() { return mat2x3(mat2(la::identity), vec2(0.0)); } + /** * Sine function where multiples of 90 degrees come out exact. * * @param x Angle in degrees. */ inline double sind(double x) { - if (!std::isfinite(x)) return sin(x); + if (!la::isfinite(x)) return sin(x); if (x < 0.0) return -sind(-x); int quo; x = remquo(fabs(x), 90.0, &quo); switch (quo % 4) { case 0: - return sin(glm::radians(x)); + return sin(radians(x)); case 1: - return cos(glm::radians(x)); + return cos(radians(x)); case 2: - return -sin(glm::radians(x)); + return -sin(radians(x)); case 3: - return -cos(glm::radians(x)); + return -cos(radians(x)); } return 0.0; } @@ -119,66 +131,65 @@ struct Box { /** * Default constructor is an infinite box that contains all space. */ - Box() {} + constexpr Box() {} /** * Creates a box that contains the two given points. */ - Box(const vec3 p1, const vec3 p2) { - min = glm::min(p1, p2); - max = glm::max(p1, p2); + constexpr Box(const vec3 p1, const vec3 p2) { + min = la::min(p1, p2); + max = la::max(p1, p2); } /** * Returns the dimensions of the Box. */ - vec3 Size() const { return max - min; } + constexpr vec3 Size() const { return max - min; } /** * Returns the center point of the Box. */ - vec3 Center() const { return 0.5 * (max + min); } + constexpr vec3 Center() const { return 0.5 * (max + min); } /** * Returns the absolute-largest coordinate value of any contained * point. */ - double Scale() const { - vec3 absMax = glm::max(glm::abs(min), glm::abs(max)); - return glm::max(absMax.x, glm::max(absMax.y, absMax.z)); + constexpr double Scale() const { + vec3 absMax = la::max(la::abs(min), la::abs(max)); + return la::max(absMax.x, la::max(absMax.y, absMax.z)); } /** * Does this box contain (includes equal) the given point? */ - bool Contains(const vec3& p) const { - return glm::all(glm::greaterThanEqual(p, min)) && - glm::all(glm::greaterThanEqual(max, p)); + constexpr bool Contains(const vec3& p) const { + return la::all(la::gequal(p, min)) && la::all(la::gequal(max, p)); } /** * Does this box contain (includes equal) the given box? */ - bool Contains(const Box& box) const { - return glm::all(glm::greaterThanEqual(box.min, min)) && - glm::all(glm::greaterThanEqual(max, box.max)); + constexpr bool Contains(const Box& box) const { + return la::all(la::gequal(box.min, min)) && + la::all(la::gequal(max, box.max)); } /** * Expand this box to include the given point. */ void Union(const vec3 p) { - min = glm::min(min, p); - max = glm::max(max, p); + min = la::min(min, p); + max = la::max(max, p); } /** * Expand this box to include the given box. */ - Box Union(const Box& box) const { + constexpr Box Union(const Box& box) const { Box out; - out.min = glm::min(min, box.min); - out.max = glm::max(max, box.max); + out.min = la::min(min, box.min); + out.max = la::max(max, box.max); return out; } @@ -189,19 +200,19 @@ struct Box { * multiples of 90 degrees), or else the resulting bounding box will no longer * bound properly. */ - Box Transform(const mat4x3& transform) const { + constexpr Box Transform(const mat3x4& transform) const { Box out; vec3 minT = transform * vec4(min, 1.0); vec3 maxT = transform * vec4(max, 1.0); - out.min = glm::min(minT, maxT); - out.max = glm::max(minT, maxT); + out.min = la::min(minT, maxT); + out.max = la::max(minT, maxT); return out; } /** * Shift this box by the given vector. */ - Box operator+(vec3 shift) const { + constexpr Box operator+(vec3 shift) const { Box out; out.min = min + shift; out.max = max + shift; @@ -220,7 +231,7 @@ struct Box { /** * Scale this box by the given vector. */ - Box operator*(vec3 scale) const { + constexpr Box operator*(vec3 scale) const { Box out; out.min = min * scale; out.max = max * scale; @@ -239,7 +250,7 @@ struct Box { /** * Does this box overlap the one given (including equality)? */ - inline bool DoesOverlap(const Box& box) const { + constexpr bool DoesOverlap(const Box& box) const { return min.x <= box.max.x && min.y <= box.max.y && min.z <= box.max.z && max.x >= box.min.x && max.y >= box.min.y && max.z >= box.min.z; } @@ -248,15 +259,15 @@ struct Box { * Does the given point project within the XY extent of this box * (including equality)? */ - inline bool DoesOverlap(vec3 p) const { // projected in z + constexpr bool DoesOverlap(vec3 p) const { // projected in z return p.x <= max.x && p.x >= min.x && p.y <= max.y && p.y >= min.y; } /** * Does this box have finite bounds? */ - bool IsFinite() const { - return glm::all(glm::isfinite(min)) && glm::all(glm::isfinite(max)); + constexpr bool IsFinite() const { + return la::all(la::isfinite(min)) && la::all(la::isfinite(max)); } }; @@ -270,14 +281,14 @@ struct Rect { /** * Default constructor is an empty rectangle.. */ - Rect() {} + constexpr Rect() {} /** * Create a rectangle that contains the two given points. */ - Rect(const vec2 a, const vec2 b) { - min = glm::min(a, b); - max = glm::max(a, b); + constexpr Rect(const vec2 a, const vec2 b) { + min = la::min(a, b); + max = la::max(a, b); } /** @name Information @@ -288,12 +299,12 @@ struct Rect { /** * Return the dimensions of the rectangle. */ - vec2 Size() const { return max - min; } + constexpr vec2 Size() const { return max - min; } /** * Return the area of the rectangle. */ - double Area() const { + constexpr double Area() const { auto sz = Size(); return sz.x * sz.y; } @@ -302,36 +313,35 @@ struct Rect { * Returns the absolute-largest coordinate value of any contained * point. */ - double Scale() const { - vec2 absMax = glm::max(glm::abs(min), glm::abs(max)); - return glm::max(absMax.x, absMax.y); + constexpr double Scale() const { + vec2 absMax = la::max(la::abs(min), la::abs(max)); + return la::max(absMax.x, absMax.y); } /** * Returns the center point of the rectangle. */ - vec2 Center() const { return 0.5 * (max + min); } + constexpr vec2 Center() const { return 0.5 * (max + min); } /** * Does this rectangle contain (includes on border) the given point? */ - bool Contains(const vec2& p) const { - return glm::all(glm::greaterThanEqual(p, min)) && - glm::all(glm::greaterThanEqual(max, p)); + constexpr bool Contains(const vec2& p) const { + return la::all(la::gequal(p, min)) && la::all(la::gequal(max, p)); } /** * Does this rectangle contain (includes equal) the given rectangle? */ - bool Contains(const Rect& rect) const { - return glm::all(glm::greaterThanEqual(rect.min, min)) && - glm::all(glm::greaterThanEqual(max, rect.max)); + constexpr bool Contains(const Rect& rect) const { + return la::all(la::gequal(rect.min, min)) && + la::all(la::gequal(max, rect.max)); } /** * Does this rectangle overlap the one given (including equality)? */ - bool DoesOverlap(const Rect& rect) const { + constexpr bool DoesOverlap(const Rect& rect) const { return min.x <= rect.max.x && min.y <= rect.max.y && max.x >= rect.min.x && max.y >= rect.min.y; } @@ -339,13 +349,13 @@ struct Rect { /** * Is the rectangle empty (containing no space)? */ - bool IsEmpty() const { return max.y <= min.y || max.x <= min.x; }; + constexpr bool IsEmpty() const { return max.y <= min.y || max.x <= min.x; }; /** * Does this recangle have finite bounds? */ - bool IsFinite() const { - return glm::all(glm::isfinite(min)) && glm::all(glm::isfinite(max)); + constexpr bool IsFinite() const { + return la::all(la::isfinite(min)) && la::all(la::isfinite(max)); } ///@} @@ -358,24 +368,24 @@ struct Rect { * Expand this rectangle (in place) to include the given point. */ void Union(const vec2 p) { - min = glm::min(min, p); - max = glm::max(max, p); + min = la::min(min, p); + max = la::max(max, p); } /** * Expand this rectangle to include the given Rect. */ - Rect Union(const Rect& rect) const { + constexpr Rect Union(const Rect& rect) const { Rect out; - out.min = glm::min(min, rect.min); - out.max = glm::max(max, rect.max); + out.min = la::min(min, rect.min); + out.max = la::max(max, rect.max); return out; } /** * Shift this rectangle by the given vector. */ - Rect operator+(const vec2 shift) const { + constexpr Rect operator+(const vec2 shift) const { Rect out; out.min = min + shift; out.max = max + shift; @@ -394,7 +404,7 @@ struct Rect { /** * Scale this rectangle by the given vector. */ - Rect operator*(const vec2 scale) const { + constexpr Rect operator*(const vec2 scale) const { Rect out; out.min = min * scale; out.max = max * scale; @@ -417,7 +427,7 @@ struct Rect { * multiples of 90 degrees), or else the resulting rectangle will no longer * bound properly. */ - Rect Transform(const mat3x2& m) const { + constexpr Rect Transform(const mat2x3& m) const { Rect rect; rect.min = m * vec3(min, 1); rect.max = m * vec3(max, 1); @@ -507,7 +517,7 @@ class Quality { static int GetCircularSegments(double radius) { if (circularSegments_ > 0) return circularSegments_; int nSegA = 360.0 / circularAngle_; - int nSegL = 2.0 * radius * glm::pi() / circularEdgeLength_; + int nSegL = 2.0 * radius * kPi / circularEdgeLength_; int nSeg = fmin(nSegA, nSegL) + 3; nSeg -= nSeg % 4; return std::max(nSeg, 3); @@ -568,5 +578,3 @@ struct ExecutionParams { }; } // namespace manifold - -#undef HOST_DEVICE diff --git a/include/manifold/cross_section.h b/include/manifold/cross_section.h index bf016692a..dab68146c 100644 --- a/include/manifold/cross_section.h +++ b/include/manifold/cross_section.h @@ -94,7 +94,7 @@ class CrossSection { CrossSection Rotate(double degrees) const; CrossSection Scale(const vec2 s) const; CrossSection Mirror(const vec2 ax) const; - CrossSection Transform(const mat3x2& m) const; + CrossSection Transform(const mat2x3& m) const; CrossSection Warp(std::function warpFunc) const; CrossSection WarpBatch(std::function)> warpFunc) const; CrossSection Simplify(double epsilon = 1e-6) const; @@ -165,7 +165,7 @@ class CrossSection { private: mutable std::shared_ptr paths_; - mutable mat3x2 transform_ = mat3x2(1.0); + mutable mat2x3 transform_ = Identity2x3(); CrossSection(std::shared_ptr paths); std::shared_ptr GetPaths() const; }; diff --git a/include/manifold/iters.h b/include/manifold/iters.h index e910acd89..8b5110ae1 100644 --- a/include/manifold/iters.h +++ b/include/manifold/iters.h @@ -34,7 +34,7 @@ struct TransformIterator { using iterator_category = typename std::iterator_traits< std::remove_const_t>::iterator_category; - TransformIterator(Iter iter, F f) : iter(iter), f(f) {} + constexpr TransformIterator(Iter iter, F f) : iter(iter), f(f) {} TransformIterator& operator=(const TransformIterator& other) { if (this == &other) return *this; @@ -43,9 +43,9 @@ struct TransformIterator { return *this; } - reference operator*() const { return f(*iter); } + constexpr reference operator*() const { return f(*iter); } - reference operator[](size_t i) const { return f(iter[i]); } + constexpr reference operator[](size_t i) const { return f(iter[i]); } // prefix increment TransformIterator& operator++() { @@ -73,7 +73,7 @@ struct TransformIterator { return old; } - TransformIterator operator+(size_t n) const { + constexpr TransformIterator operator+(size_t n) const { return TransformIterator(iter + n, f); } @@ -82,7 +82,7 @@ struct TransformIterator { return *this; } - TransformIterator operator-(size_t n) const { + constexpr TransformIterator operator-(size_t n) const { return TransformIterator(iter - n, f); } @@ -91,19 +91,23 @@ struct TransformIterator { return *this; } - bool operator==(TransformIterator other) const { return iter == other.iter; } + constexpr bool operator==(TransformIterator other) const { + return iter == other.iter; + } - bool operator!=(TransformIterator other) const { + constexpr bool operator!=(TransformIterator other) const { return !(iter == other.iter); } - bool operator<(TransformIterator other) const { return iter < other.iter; } + constexpr bool operator<(TransformIterator other) const { + return iter < other.iter; + } - difference_type operator-(TransformIterator other) const { + constexpr difference_type operator-(TransformIterator other) const { return iter - other.iter; } - operator TransformIterator() const { + constexpr operator TransformIterator() const { return TransformIterator(f, iter); } }; @@ -122,8 +126,8 @@ struct CountingIterator { constexpr CountingIterator(T counter) : counter(counter) {} - value_type operator*() const { return counter; } - value_type operator[](T i) const { return counter + i; } + constexpr value_type operator*() const { return counter; } + constexpr value_type operator[](T i) const { return counter + i; } // prefix increment CountingIterator& operator++() { @@ -151,7 +155,7 @@ struct CountingIterator { return old; } - CountingIterator operator+(T n) const { + constexpr CountingIterator operator+(T n) const { return CountingIterator(counter + n); } @@ -160,7 +164,7 @@ struct CountingIterator { return *this; } - CountingIterator operator-(T n) const { + constexpr CountingIterator operator-(T n) const { return CountingIterator(counter - n); } @@ -169,23 +173,24 @@ struct CountingIterator { return *this; } - friend bool operator==(CountingIterator a, CountingIterator b) { + constexpr friend bool operator==(CountingIterator a, CountingIterator b) { return a.counter == b.counter; } - friend bool operator!=(CountingIterator a, CountingIterator b) { + constexpr friend bool operator!=(CountingIterator a, CountingIterator b) { return a.counter != b.counter; } - friend bool operator<(CountingIterator a, CountingIterator b) { + constexpr friend bool operator<(CountingIterator a, CountingIterator b) { return a.counter < b.counter; } - friend difference_type operator-(CountingIterator a, CountingIterator b) { + constexpr friend difference_type operator-(CountingIterator a, + CountingIterator b) { return a.counter - b.counter; } - operator CountingIterator() const { + constexpr operator CountingIterator() const { return CountingIterator(counter); } }; @@ -214,15 +219,16 @@ struct StridedRange { using iterator_category = typename std::iterator_traits< std::remove_const_t>::iterator_category; - StridedRangeIter(Iter iter, int stride) : iter(iter), stride(stride) {} + constexpr StridedRangeIter(Iter iter, int stride) + : iter(iter), stride(stride) {} - reference operator*() { return *iter; } + constexpr reference operator*() { return *iter; } - std::add_const_t operator*() const { return *iter; } + constexpr std::add_const_t operator*() const { return *iter; } - reference operator[](size_t i) { return iter[i * stride]; } + constexpr reference operator[](size_t i) { return iter[i * stride]; } - std::add_const_t operator[](size_t i) const { + constexpr std::add_const_t operator[](size_t i) const { return iter[i * stride]; } @@ -252,7 +258,7 @@ struct StridedRange { return old; } - StridedRangeIter operator+(size_t n) const { + constexpr StridedRangeIter operator+(size_t n) const { return StridedRangeIter(iter + n * stride, stride); } @@ -261,7 +267,7 @@ struct StridedRange { return *this; } - StridedRangeIter operator-(size_t n) const { + constexpr StridedRangeIter operator-(size_t n) const { return StridedRangeIter(iter - n * stride, stride); } @@ -270,19 +276,20 @@ struct StridedRange { return *this; } - friend bool operator==(StridedRangeIter a, StridedRangeIter b) { + constexpr friend bool operator==(StridedRangeIter a, StridedRangeIter b) { return a.iter == b.iter; } - friend bool operator!=(StridedRangeIter a, StridedRangeIter b) { + constexpr friend bool operator!=(StridedRangeIter a, StridedRangeIter b) { return !(a.iter == b.iter); } - friend bool operator<(StridedRangeIter a, StridedRangeIter b) { + constexpr friend bool operator<(StridedRangeIter a, StridedRangeIter b) { return a.iter < b.iter; } - friend difference_type operator-(StridedRangeIter a, StridedRangeIter b) { + constexpr friend difference_type operator-(StridedRangeIter a, + StridedRangeIter b) { // note that this is not well-defined if a.stride != b.stride... return (a.iter - b.iter) / a.stride; } @@ -291,12 +298,14 @@ struct StridedRange { const size_t stride; public: - StridedRange(Iter start, Iter end, size_t stride) + constexpr StridedRange(Iter start, Iter end, size_t stride) : _start(start), _end(end), stride(stride) {} - StridedRangeIter begin() const { return StridedRangeIter(_start, stride); } + constexpr StridedRangeIter begin() const { + return StridedRangeIter(_start, stride); + } - StridedRangeIter end() const { + constexpr StridedRangeIter end() const { return StridedRangeIter(_start, stride) + ((std::distance(_start, _end) + (stride - 1)) / stride); } diff --git a/include/manifold/linalg.h b/include/manifold/linalg.h new file mode 100644 index 000000000..a26c5eb97 --- /dev/null +++ b/include/manifold/linalg.h @@ -0,0 +1,2095 @@ +// linalg.h - 2.2 - Single-header public domain linear algebra library +// +// The intent of this library is to provide the bulk of the functionality +// you need to write programs that frequently use small, fixed-size vectors +// and matrices, in domains such as computational geometry or computer +// graphics. It strives for terse, readable source code. +// +// The original author of this software is Sterling Orsten, and its permanent +// home is . If you find this software +// useful, an acknowledgement in your source text and/or product documentation +// is appreciated, but not required. +// +// The author acknowledges significant insights and contributions by: +// Stan Melax +// Dimitri Diakopoulos +// +// Some features are deprecated. Define LINALG_FORWARD_COMPATIBLE to remove +// them. + +// This is free and unencumbered software released into the public domain. +// +// Anyone is free to copy, modify, publish, use, compile, sell, or +// distribute this software, either in source code form or as a compiled +// binary, for any purpose, commercial or non-commercial, and by any +// means. +// +// In jurisdictions that recognize copyright laws, the author or authors +// of this software dedicate any and all copyright interest in the +// software to the public domain. We make this dedication for the benefit +// of the public at large and to the detriment of our heirs and +// successors. We intend this dedication to be an overt act of +// relinquishment in perpetuity of all present and future rights to this +// software under copyright law. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +// IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR +// OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, +// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// +// For more information, please refer to + +#pragma once +#ifndef LINALG_H +#define LINALG_H + +#include // For std::array +#include // For various unary math functions, such as std::sqrt +#include // For implementing namespace linalg::aliases +#include // To resolve std::abs ambiguity on clang +#include // For std::hash declaration +#include // For forward definitions of std::ostream +#include // For std::enable_if, std::is_same, std::declval + +// In Visual Studio 2015, `constexpr` applied to a member function implies +// `const`, which causes ambiguous overload resolution +#if _MSC_VER <= 1900 +#define LINALG_CONSTEXPR14 +#else +#define LINALG_CONSTEXPR14 constexpr +#endif + +namespace linalg { +// Small, fixed-length vector type, consisting of exactly M elements of type T, +// and presumed to be a column-vector unless otherwise noted. +template +struct vec; + +// Small, fixed-size matrix type, consisting of exactly M rows and N columns of +// type T, stored in column-major order. +template +struct mat; + +// Specialize converter with a function application operator that converts +// type U to type T to enable implicit conversions +template +struct converter {}; +namespace detail { +template +using conv_t = typename std::enable_if::value, + decltype(converter{}( + std::declval()))>::type; + +// Trait for retrieving scalar type of any linear algebra object +template +struct scalar_type {}; +template +struct scalar_type> { + using type = T; +}; +template +struct scalar_type> { + using type = T; +}; + +// Type returned by the compare(...) function which supports all six comparison +// operators against 0 +template +struct ord { + T a, b; +}; +template +constexpr bool operator==(const ord &o, std::nullptr_t) { + return o.a == o.b; +} +template +constexpr bool operator!=(const ord &o, std::nullptr_t) { + return !(o.a == o.b); +} +template +constexpr bool operator<(const ord &o, std::nullptr_t) { + return o.a < o.b; +} +template +constexpr bool operator>(const ord &o, std::nullptr_t) { + return o.b < o.a; +} +template +constexpr bool operator<=(const ord &o, std::nullptr_t) { + return !(o.b < o.a); +} +template +constexpr bool operator>=(const ord &o, std::nullptr_t) { + return !(o.a < o.b); +} + +// Patterns which can be used with the compare(...) function +template +struct any_compare {}; +template +struct any_compare, vec> { + using type = ord; + constexpr ord operator()(const vec &a, const vec &b) const { + return ord{a.x, b.x}; + } +}; +template +struct any_compare, vec> { + using type = ord; + constexpr ord operator()(const vec &a, const vec &b) const { + return !(a.x == b.x) ? ord{a.x, b.x} : ord{a.y, b.y}; + } +}; +template +struct any_compare, vec> { + using type = ord; + constexpr ord operator()(const vec &a, const vec &b) const { + return !(a.x == b.x) ? ord{a.x, b.x} + : !(a.y == b.y) ? ord{a.y, b.y} + : ord{a.z, b.z}; + } +}; +template +struct any_compare, vec> { + using type = ord; + constexpr ord operator()(const vec &a, const vec &b) const { + return !(a.x == b.x) ? ord{a.x, b.x} + : !(a.y == b.y) ? ord{a.y, b.y} + : !(a.z == b.z) ? ord{a.z, b.z} + : ord{a.w, b.w}; + } +}; +template +struct any_compare, mat> { + using type = ord; + constexpr ord operator()(const mat &a, + const mat &b) const { + return compare(a.x, b.x); + } +}; +template +struct any_compare, mat> { + using type = ord; + constexpr ord operator()(const mat &a, + const mat &b) const { + return a.x != b.x ? compare(a.x, b.x) : compare(a.y, b.y); + } +}; +template +struct any_compare, mat> { + using type = ord; + constexpr ord operator()(const mat &a, + const mat &b) const { + return a.x != b.x ? compare(a.x, b.x) + : a.y != b.y ? compare(a.y, b.y) + : compare(a.z, b.z); + } +}; +template +struct any_compare, mat> { + using type = ord; + constexpr ord operator()(const mat &a, + const mat &b) const { + return a.x != b.x ? compare(a.x, b.x) + : a.y != b.y ? compare(a.y, b.y) + : a.z != b.z ? compare(a.z, b.z) + : compare(a.w, b.w); + } +}; + +// Helper for compile-time index-based access to members of vector and matrix +// types +template +struct getter; +template <> +struct getter<0> { + template + constexpr auto operator()(A &a) const -> decltype(a.x) { + return a.x; + } +}; +template <> +struct getter<1> { + template + constexpr auto operator()(A &a) const -> decltype(a.y) { + return a.y; + } +}; +template <> +struct getter<2> { + template + constexpr auto operator()(A &a) const -> decltype(a.z) { + return a.z; + } +}; +template <> +struct getter<3> { + template + constexpr auto operator()(A &a) const -> decltype(a.w) { + return a.w; + } +}; + +// Stand-in for std::integer_sequence/std::make_integer_sequence +template +struct seq {}; +template +struct make_seq_impl; +template +struct make_seq_impl { + using type = seq<>; +}; +template +struct make_seq_impl { + using type = seq; +}; +template +struct make_seq_impl { + using type = seq; +}; +template +struct make_seq_impl { + using type = seq; +}; +template +struct make_seq_impl { + using type = seq; +}; +template +using make_seq = typename make_seq_impl::type; +template +vec constexpr swizzle(const vec &v, seq i) { + return {getter{}(v)...}; +} +template +mat constexpr swizzle(const mat &m, + seq i, seq j) { + return {swizzle(getter{}(m), i)...}; +} + +// SFINAE helpers to determine result of function application +template +using ret_t = decltype(std::declval()(std::declval()...)); + +// SFINAE helper which is defined if all provided types are scalars +struct empty {}; +template +struct scalars; +template <> +struct scalars<> { + using type = void; +}; +template +struct scalars : std::conditional::value, + scalars, empty>::type {}; +template +using scalars_t = typename scalars::type; + +// Helpers which indicate how apply(F, ...) should be called for various +// arguments +template +struct apply {}; // Patterns which contain only vectors or scalars +template +struct apply, vec> { + using type = vec, M>; + enum { size = M, mm = 0 }; + template + static constexpr type impl(seq, F f, const vec &a) { + return {f(getter{}(a))...}; + } +}; +template +struct apply, vec, vec> { + using type = vec, M>; + enum { size = M, mm = 0 }; + template + static constexpr type impl(seq, F f, const vec &a, + const vec &b) { + return {f(getter{}(a), getter{}(b))...}; + } +}; +template +struct apply, vec, B> { + using type = vec, M>; + enum { size = M, mm = 0 }; + template + static constexpr type impl(seq, F f, const vec &a, B b) { + return {f(getter{}(a), b)...}; + } +}; +template +struct apply, A, vec> { + using type = vec, M>; + enum { size = M, mm = 0 }; + template + static constexpr type impl(seq, F f, A a, const vec &b) { + return {f(a, getter{}(b))...}; + } +}; +template +struct apply, vec, vec, vec> { + using type = vec, M>; + enum { size = M, mm = 0 }; + template + static constexpr type impl(seq, F f, const vec &a, + const vec &b, const vec &c) { + return {f(getter{}(a), getter{}(b), getter{}(c))...}; + } +}; +template +struct apply, vec, vec, C> { + using type = vec, M>; + enum { size = M, mm = 0 }; + template + static constexpr type impl(seq, F f, const vec &a, + const vec &b, C c) { + return {f(getter{}(a), getter{}(b), c)...}; + } +}; +template +struct apply, vec, B, vec> { + using type = vec, M>; + enum { size = M, mm = 0 }; + template + static constexpr type impl(seq, F f, const vec &a, B b, + const vec &c) { + return {f(getter{}(a), b, getter{}(c))...}; + } +}; +template +struct apply, vec, B, C> { + using type = vec, M>; + enum { size = M, mm = 0 }; + template + static constexpr type impl(seq, F f, const vec &a, B b, C c) { + return {f(getter{}(a), b, c)...}; + } +}; +template +struct apply, A, vec, vec> { + using type = vec, M>; + enum { size = M, mm = 0 }; + template + static constexpr type impl(seq, F f, A a, const vec &b, + const vec &c) { + return {f(a, getter{}(b), getter{}(c))...}; + } +}; +template +struct apply, A, vec, C> { + using type = vec, M>; + enum { size = M, mm = 0 }; + template + static constexpr type impl(seq, F f, A a, const vec &b, C c) { + return {f(a, getter{}(b), c)...}; + } +}; +template +struct apply, A, B, vec> { + using type = vec, M>; + enum { size = M, mm = 0 }; + template + static constexpr type impl(seq, F f, A a, B b, const vec &c) { + return {f(a, b, getter{}(c))...}; + } +}; +template +struct apply, mat> { + using type = mat, M, N>; + enum { size = N, mm = 0 }; + template + static constexpr type impl(seq, F f, const mat &a) { + return {apply>::impl(make_seq<0, M>{}, f, + getter{}(a))...}; + } +}; +template +struct apply, mat, mat> { + using type = mat, M, N>; + enum { size = N, mm = 1 }; + template + static constexpr type impl(seq, F f, const mat &a, + const mat &b) { + return {apply, vec>::impl( + make_seq<0, M>{}, f, getter{}(a), getter{}(b))...}; + } +}; +template +struct apply, mat, B> { + using type = mat, M, N>; + enum { size = N, mm = 0 }; + template + static constexpr type impl(seq, F f, const mat &a, B b) { + return {apply, B>::impl(make_seq<0, M>{}, f, + getter{}(a), b)...}; + } +}; +template +struct apply, A, mat> { + using type = mat, M, N>; + enum { size = N, mm = 0 }; + template + static constexpr type impl(seq, F f, A a, const mat &b) { + return {apply>::impl(make_seq<0, M>{}, f, a, + getter{}(b))...}; + } +}; +template +struct apply, A...> { + using type = ret_t; + enum { size = 0, mm = 0 }; + static constexpr type impl(seq<>, F f, A... a) { return f(a...); } +}; + +// Function objects for selecting between alternatives +struct min { + template + constexpr auto operator()(A a, B b) const -> + typename std::remove_reference::type { + return a < b ? a : b; + } +}; +struct max { + template + constexpr auto operator()(A a, B b) const -> + typename std::remove_reference::type { + return a < b ? b : a; + } +}; +struct clamp { + template + constexpr auto operator()(A a, B b, C c) const -> + typename std::remove_reference::type { + return a < b ? b : a < c ? a : c; + } +}; +struct select { + template + constexpr auto operator()(A a, B b, C c) const -> + typename std::remove_reference::type { + return a ? b : c; + } +}; +struct lerp { + template + constexpr auto operator()(A a, B b, + C c) const -> decltype(a * (1 - c) + b * c) { + return a * (1 - c) + b * c; + } +}; + +// Function objects for applying operators +struct op_pos { + template + constexpr auto operator()(A a) const -> decltype(+a) { + return +a; + } +}; +struct op_neg { + template + constexpr auto operator()(A a) const -> decltype(-a) { + return -a; + } +}; +struct op_not { + template + constexpr auto operator()(A a) const -> decltype(!a) { + return !a; + } +}; +struct op_cmp { + template + constexpr auto operator()(A a) const -> decltype(~(a)) { + return ~a; + } +}; +struct op_mul { + template + constexpr auto operator()(A a, B b) const -> decltype(a * b) { + return a * b; + } +}; +struct op_div { + template + constexpr auto operator()(A a, B b) const -> decltype(a / b) { + return a / b; + } +}; +struct op_mod { + template + constexpr auto operator()(A a, B b) const -> decltype(a % b) { + return a % b; + } +}; +struct op_add { + template + constexpr auto operator()(A a, B b) const -> decltype(a + b) { + return a + b; + } +}; +struct op_sub { + template + constexpr auto operator()(A a, B b) const -> decltype(a - b) { + return a - b; + } +}; +struct op_lsh { + template + constexpr auto operator()(A a, B b) const -> decltype(a << b) { + return a << b; + } +}; +struct op_rsh { + template + constexpr auto operator()(A a, B b) const -> decltype(a >> b) { + return a >> b; + } +}; +struct op_lt { + template + constexpr auto operator()(A a, B b) const -> decltype(a < b) { + return a < b; + } +}; +struct op_gt { + template + constexpr auto operator()(A a, B b) const -> decltype(a > b) { + return a > b; + } +}; +struct op_le { + template + constexpr auto operator()(A a, B b) const -> decltype(a <= b) { + return a <= b; + } +}; +struct op_ge { + template + constexpr auto operator()(A a, B b) const -> decltype(a >= b) { + return a >= b; + } +}; +struct op_eq { + template + constexpr auto operator()(A a, B b) const -> decltype(a == b) { + return a == b; + } +}; +struct op_ne { + template + constexpr auto operator()(A a, B b) const -> decltype(a != b) { + return a != b; + } +}; +struct op_int { + template + constexpr auto operator()(A a, B b) const -> decltype(a & b) { + return a & b; + } +}; +struct op_xor { + template + constexpr auto operator()(A a, B b) const -> decltype(a ^ b) { + return a ^ b; + } +}; +struct op_un { + template + constexpr auto operator()(A a, B b) const -> decltype(a | b) { + return a | b; + } +}; +struct op_and { + template + constexpr auto operator()(A a, B b) const -> decltype(a && b) { + return a && b; + } +}; +struct op_or { + template + constexpr auto operator()(A a, B b) const -> decltype(a || b) { + return a || b; + } +}; + +// Function objects for applying standard library math functions +struct std_isfinite { + template + constexpr auto operator()(A a) const -> decltype(std::isfinite(a)) { + return std::isfinite(a); + } +}; +struct std_abs { + template + constexpr auto operator()(A a) const -> decltype(std::abs(a)) { + return std::abs(a); + } +}; +struct std_floor { + template + constexpr auto operator()(A a) const -> decltype(std::floor(a)) { + return std::floor(a); + } +}; +struct std_ceil { + template + constexpr auto operator()(A a) const -> decltype(std::ceil(a)) { + return std::ceil(a); + } +}; +struct std_exp { + template + constexpr auto operator()(A a) const -> decltype(std::exp(a)) { + return std::exp(a); + } +}; +struct std_log { + template + constexpr auto operator()(A a) const -> decltype(std::log(a)) { + return std::log(a); + } +}; +struct std_log2 { + template + constexpr auto operator()(A a) const -> decltype(std::log2(a)) { + return std::log2(a); + } +}; +struct std_log10 { + template + constexpr auto operator()(A a) const -> decltype(std::log10(a)) { + return std::log10(a); + } +}; +struct std_sqrt { + template + constexpr auto operator()(A a) const -> decltype(std::sqrt(a)) { + return std::sqrt(a); + } +}; +struct std_sin { + template + constexpr auto operator()(A a) const -> decltype(std::sin(a)) { + return std::sin(a); + } +}; +struct std_cos { + template + constexpr auto operator()(A a) const -> decltype(std::cos(a)) { + return std::cos(a); + } +}; +struct std_tan { + template + constexpr auto operator()(A a) const -> decltype(std::tan(a)) { + return std::tan(a); + } +}; +struct std_asin { + template + constexpr auto operator()(A a) const -> decltype(std::asin(a)) { + return std::asin(a); + } +}; +struct std_acos { + template + constexpr auto operator()(A a) const -> decltype(std::acos(a)) { + return std::acos(a); + } +}; +struct std_atan { + template + constexpr auto operator()(A a) const -> decltype(std::atan(a)) { + return std::atan(a); + } +}; +struct std_sinh { + template + constexpr auto operator()(A a) const -> decltype(std::sinh(a)) { + return std::sinh(a); + } +}; +struct std_cosh { + template + constexpr auto operator()(A a) const -> decltype(std::cosh(a)) { + return std::cosh(a); + } +}; +struct std_tanh { + template + constexpr auto operator()(A a) const -> decltype(std::tanh(a)) { + return std::tanh(a); + } +}; +struct std_round { + template + constexpr auto operator()(A a) const -> decltype(std::round(a)) { + return std::round(a); + } +}; +struct std_fmod { + template + constexpr auto operator()(A a, B b) const -> decltype(std::fmod(a, b)) { + return std::fmod(a, b); + } +}; +struct std_pow { + template + constexpr auto operator()(A a, B b) const -> decltype(std::pow(a, b)) { + return std::pow(a, b); + } +}; +struct std_atan2 { + template + constexpr auto operator()(A a, B b) const -> decltype(std::atan2(a, b)) { + return std::atan2(a, b); + } +}; +struct std_copysign { + template + constexpr auto operator()(A a, B b) const -> decltype(std::copysign(a, b)) { + return std::copysign(a, b); + } +}; +} // namespace detail + +// Small, fixed-length vector type, consisting of exactly M elements of type T, +// and presumed to be a column-vector unless otherwise noted +template +struct vec { + T x; + constexpr vec() : x() {} + constexpr vec(const T &x_) : x(x_) {} + // NOTE: vec does NOT have a constructor from pointer, this can conflict + // with initializing its single element from zero + template + constexpr explicit vec(const vec &v) : vec(static_cast(v.x)) {} + constexpr const T &operator[](int i) const { return x; } + LINALG_CONSTEXPR14 T &operator[](int i) { return x; } + + template > + constexpr vec(const U &u) : vec(converter{}(u)) {} + template > + constexpr operator U() const { + return converter{}(*this); + } +}; +template +struct vec { + T x, y; + constexpr vec() : x(), y() {} + constexpr vec(const T &x_, const T &y_) : x(x_), y(y_) {} + constexpr explicit vec(const T &s) : vec(s, s) {} + constexpr explicit vec(const T *p) : vec(p[0], p[1]) {} + template + constexpr explicit vec(const vec &v) + : vec(static_cast(v.x), static_cast(v.y)) { + static_assert( + N >= 2, + "You must give extra arguments if your input vector is shorter."); + } + constexpr const T &operator[](int i) const { return i == 0 ? x : y; } + LINALG_CONSTEXPR14 T &operator[](int i) { return i == 0 ? x : y; } + + template > + constexpr vec(const U &u) : vec(converter{}(u)) {} + template > + constexpr operator U() const { + return converter{}(*this); + } +}; +template +struct vec { + T x, y, z; + constexpr vec() : x(), y(), z() {} + constexpr vec(const T &x_, const T &y_, const T &z_) : x(x_), y(y_), z(z_) {} + constexpr vec(const vec &xy, const T &z_) : vec(xy.x, xy.y, z_) {} + constexpr explicit vec(const T &s) : vec(s, s, s) {} + constexpr explicit vec(const T *p) : vec(p[0], p[1], p[2]) {} + template + constexpr explicit vec(const vec &v) + : vec(static_cast(v.x), static_cast(v.y), static_cast(v.z)) { + static_assert( + N >= 3, + "You must give extra arguments if your input vector is shorter."); + } + constexpr const T &operator[](int i) const { + return i == 0 ? x : i == 1 ? y : z; + } + LINALG_CONSTEXPR14 T &operator[](int i) { + return i == 0 ? x : i == 1 ? y : z; + } + constexpr const vec &xy() const { + return *reinterpret_cast *>(this); + } + vec &xy() { return *reinterpret_cast *>(this); } + + template > + constexpr vec(const U &u) : vec(converter{}(u)) {} + template > + constexpr operator U() const { + return converter{}(*this); + } +}; +template +struct vec { + T x, y, z, w; + constexpr vec() : x(), y(), z(), w() {} + constexpr vec(const T &x_, const T &y_, const T &z_, const T &w_) + : x(x_), y(y_), z(z_), w(w_) {} + constexpr vec(const vec &xy, const T &z_, const T &w_) + : vec(xy.x, xy.y, z_, w_) {} + constexpr vec(const vec &xyz, const T &w_) + : vec(xyz.x, xyz.y, xyz.z, w_) {} + constexpr explicit vec(const T &s) : vec(s, s, s, s) {} + constexpr explicit vec(const T *p) : vec(p[0], p[1], p[2], p[3]) {} + template + constexpr explicit vec(const vec &v) + : vec(static_cast(v.x), static_cast(v.y), static_cast(v.z), + static_cast(v.w)) { + static_assert( + N >= 4, + "You must give extra arguments if your input vector is shorter."); + } + constexpr const T &operator[](int i) const { + return i == 0 ? x : i == 1 ? y : i == 2 ? z : w; + } + LINALG_CONSTEXPR14 T &operator[](int i) { + return i == 0 ? x : i == 1 ? y : i == 2 ? z : w; + } + constexpr const vec &xy() const { + return *reinterpret_cast *>(this); + } + constexpr const vec &xyz() const { + return *reinterpret_cast *>(this); + } + vec &xy() { return *reinterpret_cast *>(this); } + vec &xyz() { return *reinterpret_cast *>(this); } + + template > + constexpr vec(const U &u) : vec(converter{}(u)) {} + template > + constexpr operator U() const { + return converter{}(*this); + } +}; + +// Small, fixed-size matrix type, consisting of exactly M rows and N columns of +// type T, stored in column-major order. +template +struct mat { + typedef vec V; + V x; + constexpr mat() : x() {} + constexpr mat(const V &x_) : x(x_) {} + constexpr explicit mat(const T &s) : x(s) {} + constexpr explicit mat(const T *p) : x(p + M * 0) {} + template + constexpr explicit mat(const mat &m) : mat(V(m.x)) {} + constexpr vec row(int i) const { return {x[i]}; } + constexpr const V &operator[](int j) const { return x; } + LINALG_CONSTEXPR14 V &operator[](int j) { return x; } + + template > + constexpr mat(const U &u) : mat(converter{}(u)) {} + template > + constexpr operator U() const { + return converter{}(*this); + } +}; +template +struct mat { + typedef vec V; + V x, y; + constexpr mat() : x(), y() {} + constexpr mat(const V &x_, const V &y_) : x(x_), y(y_) {} + constexpr explicit mat(const T &s) : x(s), y(s) {} + constexpr explicit mat(const T *p) : x(p + M * 0), y(p + M * 1) {} + template + constexpr explicit mat(const mat &m) : mat(V(m.x), V(m.y)) { + static_assert(P >= 2, "Input matrix dimensions must be at least as big."); + } + constexpr vec row(int i) const { return {x[i], y[i]}; } + constexpr const V &operator[](int j) const { return j == 0 ? x : y; } + LINALG_CONSTEXPR14 V &operator[](int j) { return j == 0 ? x : y; } + + template > + constexpr mat(const U &u) : mat(converter{}(u)) {} + template > + constexpr operator U() const { + return converter{}(*this); + } +}; +template +struct mat { + typedef vec V; + V x, y, z; + constexpr mat() : x(), y(), z() {} + constexpr mat(const V &x_, const V &y_, const V &z_) : x(x_), y(y_), z(z_) {} + constexpr mat(const mat &m_, const V &z_) + : x(m_.x), y(m_.y), z(z_) {} + constexpr explicit mat(const T &s) : x(s), y(s), z(s) {} + constexpr explicit mat(const T *p) + : x(p + M * 0), y(p + M * 1), z(p + M * 2) {} + template + constexpr explicit mat(const mat &m) : mat(V(m.x), V(m.y), V(m.z)) { + static_assert(P >= 3, "Input matrix dimensions must be at least as big."); + } + constexpr vec row(int i) const { return {x[i], y[i], z[i]}; } + constexpr const V &operator[](int j) const { + return j == 0 ? x : j == 1 ? y : z; + } + LINALG_CONSTEXPR14 V &operator[](int j) { + return j == 0 ? x : j == 1 ? y : z; + } + + template > + constexpr mat(const U &u) : mat(converter{}(u)) {} + template > + constexpr operator U() const { + return converter{}(*this); + } +}; +template +struct mat { + typedef vec V; + V x, y, z, w; + constexpr mat() : x(), y(), z(), w() {} + constexpr mat(const V &x_, const V &y_, const V &z_, const V &w_) + : x(x_), y(y_), z(z_), w(w_) {} + constexpr mat(const mat &m_, const V &w_) + : x(m_.x), y(m_.y), z(m_.z), w(w_) {} + constexpr explicit mat(const T &s) : x(s), y(s), z(s), w(s) {} + constexpr explicit mat(const T *p) + : x(p + M * 0), y(p + M * 1), z(p + M * 2), w(p + M * 3) {} + template + constexpr explicit mat(const mat &m) + : mat(V(m.x), V(m.y), V(m.z), V(m.w)) { + static_assert(P >= 4, "Input matrix dimensions must be at least as big."); + } + + constexpr vec row(int i) const { return {x[i], y[i], z[i], w[i]}; } + constexpr const V &operator[](int j) const { + return j == 0 ? x : j == 1 ? y : j == 2 ? z : w; + } + LINALG_CONSTEXPR14 V &operator[](int j) { + return j == 0 ? x : j == 1 ? y : j == 2 ? z : w; + } + + template > + constexpr mat(const U &u) : mat(converter{}(u)) {} + template > + constexpr operator U() const { + return converter{}(*this); + } +}; + +// Define a type which will convert to the multiplicative identity of any square +// matrix +struct identity_t { + constexpr explicit identity_t(int) {} +}; +template +struct converter, identity_t> { + constexpr mat operator()(identity_t) const { return {vec{1}}; } +}; +template +struct converter, identity_t> { + constexpr mat operator()(identity_t) const { + return {{1, 0}, {0, 1}}; + } +}; +template +struct converter, identity_t> { + constexpr mat operator()(identity_t) const { + return {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; + } +}; +template +struct converter, identity_t> { + constexpr mat operator()(identity_t) const { + return {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}; + } +}; +constexpr identity_t identity{1}; + +// Produce a scalar by applying f(A,B) -> A to adjacent pairs of elements from a +// vec/mat in left-to-right/column-major order (matching the associativity of +// arithmetic and logical operators) +template +constexpr A fold(F f, A a, const vec &b) { + return f(a, b.x); +} +template +constexpr A fold(F f, A a, const vec &b) { + return f(f(a, b.x), b.y); +} +template +constexpr A fold(F f, A a, const vec &b) { + return f(f(f(a, b.x), b.y), b.z); +} +template +constexpr A fold(F f, A a, const vec &b) { + return f(f(f(f(a, b.x), b.y), b.z), b.w); +} +template +constexpr A fold(F f, A a, const mat &b) { + return fold(f, a, b.x); +} +template +constexpr A fold(F f, A a, const mat &b) { + return fold(f, fold(f, a, b.x), b.y); +} +template +constexpr A fold(F f, A a, const mat &b) { + return fold(f, fold(f, fold(f, a, b.x), b.y), b.z); +} +template +constexpr A fold(F f, A a, const mat &b) { + return fold(f, fold(f, fold(f, fold(f, a, b.x), b.y), b.z), b.w); +} + +// Type aliases for the result of calling apply(...) with various arguments, can +// be used with return type SFINAE to constrian overload sets +template +using apply_t = typename detail::apply::type; +template +using mm_apply_t = typename std::enable_if::mm, + apply_t>::type; +template +using no_mm_apply_t = typename std::enable_if::mm, + apply_t>::type; +template +using scalar_t = + typename detail::scalar_type::type; // Underlying scalar type when + // performing elementwise operations + +// apply(f,...) applies the provided function in an elementwise fashion to its +// arguments, producing an object of the same dimensions +template +constexpr apply_t apply(F func, const A &...args) { + return detail::apply::impl( + detail::make_seq<0, detail::apply::size>{}, func, args...); +} + +// map(a,f) is equivalent to apply(f,a) +template +constexpr apply_t map(const A &a, F func) { + return apply(func, a); +} + +// zip(a,b,f) is equivalent to apply(f,a,b) +template +constexpr apply_t zip(const A &a, const B &b, F func) { + return apply(func, a, b); +} + +// Relational operators are defined to compare the elements of two vectors or +// matrices lexicographically, in column-major order +template +constexpr typename detail::any_compare::type compare(const A &a, + const B &b) { + return detail::any_compare()(a, b); +} +template +constexpr auto operator==(const A &a, + const B &b) -> decltype(compare(a, b) == 0) { + return compare(a, b) == 0; +} +template +constexpr auto operator!=(const A &a, + const B &b) -> decltype(compare(a, b) != 0) { + return compare(a, b) != 0; +} +template +constexpr auto operator<(const A &a, + const B &b) -> decltype(compare(a, b) < 0) { + return compare(a, b) < 0; +} +template +constexpr auto operator>(const A &a, + const B &b) -> decltype(compare(a, b) > 0) { + return compare(a, b) > 0; +} +template +constexpr auto operator<=(const A &a, + const B &b) -> decltype(compare(a, b) <= 0) { + return compare(a, b) <= 0; +} +template +constexpr auto operator>=(const A &a, + const B &b) -> decltype(compare(a, b) >= 0) { + return compare(a, b) >= 0; +} + +// Functions for coalescing scalar values +template +constexpr bool any(const A &a) { + return fold(detail::op_or{}, false, a); +} +template +constexpr bool all(const A &a) { + return fold(detail::op_and{}, true, a); +} +template +constexpr scalar_t sum(const A &a) { + return fold(detail::op_add{}, scalar_t(0), a); +} +template +constexpr scalar_t product(const A &a) { + return fold(detail::op_mul{}, scalar_t(1), a); +} +template +constexpr scalar_t minelem(const A &a) { + return fold(detail::min{}, a.x, a); +} +template +constexpr scalar_t maxelem(const A &a) { + return fold(detail::max{}, a.x, a); +} +template +int argmin(const vec &a) { + int j = 0; + for (int i = 1; i < M; ++i) + if (a[i] < a[j]) j = i; + return j; +} +template +int argmax(const vec &a) { + int j = 0; + for (int i = 1; i < M; ++i) + if (a[i] > a[j]) j = i; + return j; +} + +// Unary operators are defined component-wise for linalg types +template +constexpr apply_t operator+(const A &a) { + return apply(detail::op_pos{}, a); +} +template +constexpr apply_t operator-(const A &a) { + return apply(detail::op_neg{}, a); +} +template +constexpr apply_t operator~(const A &a) { + return apply(detail::op_cmp{}, a); +} +template +constexpr apply_t operator!(const A &a) { + return apply(detail::op_not{}, a); +} + +// Binary operators are defined component-wise for linalg types, EXCEPT for +// `operator *` +template +constexpr apply_t operator+(const A &a, const B &b) { + return apply(detail::op_add{}, a, b); +} +template +constexpr apply_t operator-(const A &a, const B &b) { + return apply(detail::op_sub{}, a, b); +} +template +constexpr apply_t cmul(const A &a, const B &b) { + return apply(detail::op_mul{}, a, b); +} +template +constexpr apply_t operator/(const A &a, const B &b) { + return apply(detail::op_div{}, a, b); +} +template +constexpr apply_t operator%(const A &a, const B &b) { + return apply(detail::op_mod{}, a, b); +} +template +constexpr apply_t operator|(const A &a, const B &b) { + return apply(detail::op_un{}, a, b); +} +template +constexpr apply_t operator^(const A &a, const B &b) { + return apply(detail::op_xor{}, a, b); +} +template +constexpr apply_t operator&(const A &a, const B &b) { + return apply(detail::op_int{}, a, b); +} +template +constexpr apply_t operator<<(const A &a, const B &b) { + return apply(detail::op_lsh{}, a, b); +} +template +constexpr apply_t operator>>(const A &a, const B &b) { + return apply(detail::op_rsh{}, a, b); +} + +// Binary `operator *` represents the algebraic matrix product - use cmul(a, b) +// for the Hadamard (component-wise) product. +template +constexpr auto operator*(const A &a, const B &b) { + return mul(a, b); +} + +// Binary assignment operators a $= b is always defined as though it were +// explicitly written a = a $ b +template +constexpr auto operator+=(A &a, const B &b) -> decltype(a = a + b) { + return a = a + b; +} +template +constexpr auto operator-=(A &a, const B &b) -> decltype(a = a - b) { + return a = a - b; +} +template +constexpr auto operator*=(A &a, const B &b) -> decltype(a = a * b) { + return a = a * b; +} +template +constexpr auto operator/=(A &a, const B &b) -> decltype(a = a / b) { + return a = a / b; +} +template +constexpr auto operator%=(A &a, const B &b) -> decltype(a = a % b) { + return a = a % b; +} +template +constexpr auto operator|=(A &a, const B &b) -> decltype(a = a | b) { + return a = a | b; +} +template +constexpr auto operator^=(A &a, const B &b) -> decltype(a = a ^ b) { + return a = a ^ b; +} +template +constexpr auto operator&=(A &a, const B &b) -> decltype(a = a & b) { + return a = a & b; +} +template +constexpr auto operator<<=(A &a, const B &b) -> decltype(a = a << b) { + return a = a << b; +} +template +constexpr auto operator>>=(A &a, const B &b) -> decltype(a = a >> b) { + return a = a >> b; +} + +// Swizzles and subobjects +template +constexpr vec swizzle(const vec &a) { + return {detail::getter{}(a)...}; +} +template +constexpr vec subvec(const vec &a) { + return detail::swizzle(a, detail::make_seq{}); +} +template +constexpr mat submat(const mat &a) { + return detail::swizzle(a, detail::make_seq{}, + detail::make_seq{}); +} + +// Component-wise standard library math functions +template +constexpr apply_t isfinite(const A &a) { + return apply(detail::std_isfinite{}, a); +} +template +constexpr apply_t abs(const A &a) { + return apply(detail::std_abs{}, a); +} +template +constexpr apply_t floor(const A &a) { + return apply(detail::std_floor{}, a); +} +template +constexpr apply_t ceil(const A &a) { + return apply(detail::std_ceil{}, a); +} +template +constexpr apply_t exp(const A &a) { + return apply(detail::std_exp{}, a); +} +template +constexpr apply_t log(const A &a) { + return apply(detail::std_log{}, a); +} +template +constexpr apply_t log2(const A &a) { + return apply(detail::std_log2{}, a); +} +template +constexpr apply_t log10(const A &a) { + return apply(detail::std_log10{}, a); +} +template +constexpr apply_t sqrt(const A &a) { + return apply(detail::std_sqrt{}, a); +} +template +constexpr apply_t sin(const A &a) { + return apply(detail::std_sin{}, a); +} +template +constexpr apply_t cos(const A &a) { + return apply(detail::std_cos{}, a); +} +template +constexpr apply_t tan(const A &a) { + return apply(detail::std_tan{}, a); +} +template +constexpr apply_t asin(const A &a) { + return apply(detail::std_asin{}, a); +} +template +constexpr apply_t acos(const A &a) { + return apply(detail::std_acos{}, a); +} +template +constexpr apply_t atan(const A &a) { + return apply(detail::std_atan{}, a); +} +template +constexpr apply_t sinh(const A &a) { + return apply(detail::std_sinh{}, a); +} +template +constexpr apply_t cosh(const A &a) { + return apply(detail::std_cosh{}, a); +} +template +constexpr apply_t tanh(const A &a) { + return apply(detail::std_tanh{}, a); +} +template +constexpr apply_t round(const A &a) { + return apply(detail::std_round{}, a); +} + +template +constexpr apply_t fmod(const A &a, const B &b) { + return apply(detail::std_fmod{}, a, b); +} +template +constexpr apply_t pow(const A &a, const B &b) { + return apply(detail::std_pow{}, a, b); +} +template +constexpr apply_t atan2(const A &a, const B &b) { + return apply(detail::std_atan2{}, a, b); +} +template +constexpr apply_t copysign(const A &a, const B &b) { + return apply(detail::std_copysign{}, a, b); +} + +// Component-wise relational functions on vectors +template +constexpr apply_t equal(const A &a, const B &b) { + return apply(detail::op_eq{}, a, b); +} +template +constexpr apply_t nequal(const A &a, const B &b) { + return apply(detail::op_ne{}, a, b); +} +template +constexpr apply_t less(const A &a, const B &b) { + return apply(detail::op_lt{}, a, b); +} +template +constexpr apply_t greater(const A &a, const B &b) { + return apply(detail::op_gt{}, a, b); +} +template +constexpr apply_t lequal(const A &a, const B &b) { + return apply(detail::op_le{}, a, b); +} +template +constexpr apply_t gequal(const A &a, const B &b) { + return apply(detail::op_ge{}, a, b); +} + +// Component-wise selection functions on vectors +template +constexpr apply_t min(const A &a, const B &b) { + return apply(detail::min{}, a, b); +} +template +constexpr apply_t max(const A &a, const B &b) { + return apply(detail::max{}, a, b); +} +template +constexpr apply_t clamp(const X &x, const L &l, + const H &h) { + return apply(detail::clamp{}, x, l, h); +} +template +constexpr apply_t select(const P &p, const A &a, + const B &b) { + return apply(detail::select{}, p, a, b); +} +template +constexpr apply_t lerp(const A &a, const B &b, + const T &t) { + return apply(detail::lerp{}, a, b, t); +} + +// Support for vector algebra +template +constexpr T cross(const vec &a, const vec &b) { + return a.x * b.y - a.y * b.x; +} +template +constexpr vec cross(T a, const vec &b) { + return {-a * b.y, a * b.x}; +} +template +constexpr vec cross(const vec &a, T b) { + return {a.y * b, -a.x * b}; +} +template +constexpr vec cross(const vec &a, const vec &b) { + return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x}; +} +template +constexpr T dot(const vec &a, const vec &b) { + return sum(a * b); +} +template +constexpr T length2(const vec &a) { + return dot(a, a); +} +template +T length(const vec &a) { + return std::sqrt(length2(a)); +} +template +vec normalize(const vec &a) { + return a / length(a); +} +template +constexpr T distance2(const vec &a, const vec &b) { + return length2(b - a); +} +template +T distance(const vec &a, const vec &b) { + return length(b - a); +} +template +T uangle(const vec &a, const vec &b) { + T d = dot(a, b); + return d > 1 ? 0 : std::acos(d < -1 ? -1 : d); +} +template +T angle(const vec &a, const vec &b) { + return uangle(normalize(a), normalize(b)); +} +template +vec rot(T a, const vec &v) { + const T s = std::sin(a), c = std::cos(a); + return {v.x * c - v.y * s, v.x * s + v.y * c}; +} +template +vec rotx(T a, const vec &v) { + const T s = std::sin(a), c = std::cos(a); + return {v.x, v.y * c - v.z * s, v.y * s + v.z * c}; +} +template +vec roty(T a, const vec &v) { + const T s = std::sin(a), c = std::cos(a); + return {v.x * c + v.z * s, v.y, -v.x * s + v.z * c}; +} +template +vec rotz(T a, const vec &v) { + const T s = std::sin(a), c = std::cos(a); + return {v.x * c - v.y * s, v.x * s + v.y * c, v.z}; +} +template +vec nlerp(const vec &a, const vec &b, T t) { + return normalize(lerp(a, b, t)); +} +template +vec slerp(const vec &a, const vec &b, T t) { + T th = uangle(a, b); + return th == 0 ? a + : a * (std::sin(th * (1 - t)) / std::sin(th)) + + b * (std::sin(th * t) / std::sin(th)); +} + +// Support for quaternion algebra using 4D vectors, representing xi + yj + zk + +// w +template +constexpr vec qconj(const vec &q) { + return {-q.x, -q.y, -q.z, q.w}; +} +template +vec qinv(const vec &q) { + return qconj(q) / length2(q); +} +template +vec qexp(const vec &q) { + const auto v = q.xyz(); + const auto vv = length(v); + return std::exp(q.w) * + vec{v * (vv > 0 ? std::sin(vv) / vv : 0), std::cos(vv)}; +} +template +vec qlog(const vec &q) { + const auto v = q.xyz(); + const auto vv = length(v), qq = length(q); + return {v * (vv > 0 ? std::acos(q.w / qq) / vv : 0), std::log(qq)}; +} +template +vec qpow(const vec &q, const T &p) { + const auto v = q.xyz(); + const auto vv = length(v), qq = length(q), th = std::acos(q.w / qq); + return std::pow(qq, p) * + vec{v * (vv > 0 ? std::sin(p * th) / vv : 0), std::cos(p * th)}; +} +template +constexpr vec qmul(const vec &a, const vec &b) { + return {a.x * b.w + a.w * b.x + a.y * b.z - a.z * b.y, + a.y * b.w + a.w * b.y + a.z * b.x - a.x * b.z, + a.z * b.w + a.w * b.z + a.x * b.y - a.y * b.x, + a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z}; +} +template +constexpr vec qmul(const vec &a, R... r) { + return qmul(a, qmul(r...)); +} + +// Support for 3D spatial rotations using quaternions, via qmul(qmul(q, v), +// qconj(q)) +template +constexpr vec qxdir(const vec &q) { + return {q.w * q.w + q.x * q.x - q.y * q.y - q.z * q.z, + (q.x * q.y + q.z * q.w) * 2, (q.z * q.x - q.y * q.w) * 2}; +} +template +constexpr vec qydir(const vec &q) { + return {(q.x * q.y - q.z * q.w) * 2, + q.w * q.w - q.x * q.x + q.y * q.y - q.z * q.z, + (q.y * q.z + q.x * q.w) * 2}; +} +template +constexpr vec qzdir(const vec &q) { + return {(q.z * q.x + q.y * q.w) * 2, (q.y * q.z - q.x * q.w) * 2, + q.w * q.w - q.x * q.x - q.y * q.y + q.z * q.z}; +} +template +constexpr mat qmat(const vec &q) { + return {qxdir(q), qydir(q), qzdir(q)}; +} +template +constexpr vec qrot(const vec &q, const vec &v) { + return qxdir(q) * v.x + qydir(q) * v.y + qzdir(q) * v.z; +} +template +T qangle(const vec &q) { + return std::atan2(length(q.xyz()), q.w) * 2; +} +template +vec qaxis(const vec &q) { + return normalize(q.xyz()); +} +template +vec qnlerp(const vec &a, const vec &b, T t) { + return nlerp(a, dot(a, b) < 0 ? -b : b, t); +} +template +vec qslerp(const vec &a, const vec &b, T t) { + return slerp(a, dot(a, b) < 0 ? -b : b, t); +} + +// Support for matrix algebra +template +constexpr vec mul(const vec &a, const T &b) { + return cmul(a, b); +} +template +constexpr vec mul(const T &b, const vec &a) { + return cmul(b, a); +} +template +constexpr mat mul(const mat &a, const T &b) { + return cmul(a, b); +} +template +constexpr mat mul(const T &b, const mat &a) { + return cmul(b, a); +} +template +constexpr vec mul(const vec &a, const vec &b) { + return cmul(a, b); +} +template +constexpr vec mul(const mat &a, const vec &b) { + return a.x * b.x; +} +template +constexpr vec mul(const mat &a, const vec &b) { + return a.x * b.x + a.y * b.y; +} +template +constexpr vec mul(const mat &a, const vec &b) { + return a.x * b.x + a.y * b.y + a.z * b.z; +} +template +constexpr vec mul(const mat &a, const vec &b) { + return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w; +} +template +constexpr mat mul(const mat &a, const mat &b) { + return {mul(a, b.x)}; +} +template +constexpr mat mul(const mat &a, const mat &b) { + return {mul(a, b.x), mul(a, b.y)}; +} +template +constexpr mat mul(const mat &a, const mat &b) { + return {mul(a, b.x), mul(a, b.y), mul(a, b.z)}; +} +template +constexpr mat mul(const mat &a, const mat &b) { + return {mul(a, b.x), mul(a, b.y), mul(a, b.z), mul(a, b.w)}; +} +template +constexpr vec mul(const mat &a, const mat &b, + const vec &c) { + return mul(mul(a, b), c); +} +template +constexpr mat mul(const mat &a, const mat &b, + const mat &c) { + return mul(mul(a, b), c); +} +template +constexpr vec mul(const mat &a, const mat &b, + const mat &c, const vec &d) { + return mul(mul(a, b, c), d); +} +template +constexpr mat mul(const mat &a, const mat &b, + const mat &c, const mat &d) { + return mul(mul(a, b, c), d); +} +template +constexpr mat outerprod(const vec &a, const vec &b) { + return {a * b.x}; +} +template +constexpr mat outerprod(const vec &a, const vec &b) { + return {a * b.x, a * b.y}; +} +template +constexpr mat outerprod(const vec &a, const vec &b) { + return {a * b.x, a * b.y, a * b.z}; +} +template +constexpr mat outerprod(const vec &a, const vec &b) { + return {a * b.x, a * b.y, a * b.z, a * b.w}; +} +template +constexpr vec diagonal(const mat &a) { + return {a.x.x}; +} +template +constexpr vec diagonal(const mat &a) { + return {a.x.x, a.y.y}; +} +template +constexpr vec diagonal(const mat &a) { + return {a.x.x, a.y.y, a.z.z}; +} +template +constexpr vec diagonal(const mat &a) { + return {a.x.x, a.y.y, a.z.z, a.w.w}; +} +template +constexpr T trace(const mat &a) { + return sum(diagonal(a)); +} +template +constexpr mat transpose(const mat &m) { + return {m.row(0)}; +} +template +constexpr mat transpose(const mat &m) { + return {m.row(0), m.row(1)}; +} +template +constexpr mat transpose(const mat &m) { + return {m.row(0), m.row(1), m.row(2)}; +} +template +constexpr mat transpose(const mat &m) { + return {m.row(0), m.row(1), m.row(2), m.row(3)}; +} +template +constexpr mat transpose(const vec &m) { + return transpose(mat(m)); +} +template +constexpr mat adjugate(const mat &a) { + return {vec{1}}; +} +template +constexpr mat adjugate(const mat &a) { + return {{a.y.y, -a.x.y}, {-a.y.x, a.x.x}}; +} +template +constexpr mat adjugate(const mat &a); +template +constexpr mat adjugate(const mat &a); +template +constexpr mat comatrix(const mat &a) { + return transpose(adjugate(a)); +} +template +constexpr T determinant(const mat &a) { + return a.x.x; +} +template +constexpr T determinant(const mat &a) { + return a.x.x * a.y.y - a.x.y * a.y.x; +} +template +constexpr T determinant(const mat &a) { + return a.x.x * (a.y.y * a.z.z - a.z.y * a.y.z) + + a.x.y * (a.y.z * a.z.x - a.z.z * a.y.x) + + a.x.z * (a.y.x * a.z.y - a.z.x * a.y.y); +} +template +constexpr T determinant(const mat &a); +template +constexpr mat inverse(const mat &a) { + return adjugate(a) / determinant(a); +} + +// Vectors and matrices can be used as ranges +template +T *begin(vec &a) { + return &a.x; +} +template +const T *begin(const vec &a) { + return &a.x; +} +template +T *end(vec &a) { + return begin(a) + M; +} +template +const T *end(const vec &a) { + return begin(a) + M; +} +template +vec *begin(mat &a) { + return &a.x; +} +template +const vec *begin(const mat &a) { + return &a.x; +} +template +vec *end(mat &a) { + return begin(a) + N; +} +template +const vec *end(const mat &a) { + return begin(a) + N; +} + +// Factory functions for 3D spatial transformations (will possibly be removed or +// changed in a future version) +enum fwd_axis { + neg_z, + pos_z +}; // Should projection matrices be generated assuming forward is {0,0,-1} or + // {0,0,1} +enum z_range { + neg_one_to_one, + zero_to_one +}; // Should projection matrices map z into the range of [-1,1] or [0,1]? +template +vec constexpr rotation_quat(const vec &axis, T angle) { + return {axis * std::sin(angle / 2), std::cos(angle / 2)}; +} +template +vec constexpr rotation_quat(const vec &orig, + const vec &dest) { + T cosTheta = dot(orig, dest); + if (cosTheta >= 1 - std::numeric_limits::epsilon()) { + return {0, 0, 0, 1}; + } + if (cosTheta < -1 + std::numeric_limits::epsilon()) { + vec axis = cross(vec(0, 0, 1), orig); + if (length2(axis) < std::numeric_limits::epsilon()) + axis = cross(vec(1, 0, 0), orig); + return rotation_quat(normalize(axis), + 3.14159265358979323846264338327950288); + } + vec axis = cross(orig, dest); + T s = std::sqrt((1 + cosTheta) * 2); + return {axis * (1 / s), s * 0.5}; +} +template +vec rotation_quat(const mat &m); +template +mat translation_matrix(const vec &translation) { + return {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {translation, 1}}; +} +template +mat rotation_matrix(const vec &rotation) { + return {{qxdir(rotation), 0}, + {qydir(rotation), 0}, + {qzdir(rotation), 0}, + {0, 0, 0, 1}}; +} +template +mat scaling_matrix(const vec &scaling) { + return {{scaling.x, 0, 0, 0}, + {0, scaling.y, 0, 0}, + {0, 0, scaling.z, 0}, + {0, 0, 0, 1}}; +} +template +mat pose_matrix(const vec &q, const vec &p) { + return {{qxdir(q), 0}, {qydir(q), 0}, {qzdir(q), 0}, {p, 1}}; +} +template +mat lookat_matrix(const vec &eye, const vec ¢er, + const vec &view_y_dir, fwd_axis fwd = neg_z); +template +mat frustum_matrix(T x0, T x1, T y0, T y1, T n, T f, + fwd_axis a = neg_z, z_range z = neg_one_to_one); +template +mat perspective_matrix(T fovy, T aspect, T n, T f, fwd_axis a = neg_z, + z_range z = neg_one_to_one) { + T y = n * std::tan(fovy / 2), x = y * aspect; + return frustum_matrix(-x, x, -y, y, n, f, a, z); +} + +// Provide implicit conversion between linalg::vec and std::array +template +struct converter, std::array> { + vec operator()(const std::array &a) const { return {a[0]}; } +}; +template +struct converter, std::array> { + vec operator()(const std::array &a) const { return {a[0], a[1]}; } +}; +template +struct converter, std::array> { + vec operator()(const std::array &a) const { + return {a[0], a[1], a[2]}; + } +}; +template +struct converter, std::array> { + vec operator()(const std::array &a) const { + return {a[0], a[1], a[2], a[3]}; + } +}; + +template +struct converter, vec> { + std::array operator()(const vec &a) const { return {a[0]}; } +}; +template +struct converter, vec> { + std::array operator()(const vec &a) const { return {a[0], a[1]}; } +}; +template +struct converter, vec> { + std::array operator()(const vec &a) const { + return {a[0], a[1], a[2]}; + } +}; +template +struct converter, vec> { + std::array operator()(const vec &a) const { + return {a[0], a[1], a[2], a[3]}; + } +}; +} // namespace linalg + +namespace std { +// Provide specializations for std::hash<...> with linalg types +template +struct hash> { + std::size_t operator()(const linalg::vec &v) const { + std::hash h; + return h(v.x); + } +}; +template +struct hash> { + std::size_t operator()(const linalg::vec &v) const { + std::hash h; + return h(v.x) ^ (h(v.y) << 1); + } +}; +template +struct hash> { + std::size_t operator()(const linalg::vec &v) const { + std::hash h; + return h(v.x) ^ (h(v.y) << 1) ^ (h(v.z) << 2); + } +}; +template +struct hash> { + std::size_t operator()(const linalg::vec &v) const { + std::hash h; + return h(v.x) ^ (h(v.y) << 1) ^ (h(v.z) << 2) ^ (h(v.w) << 3); + } +}; + +template +struct hash> { + std::size_t operator()(const linalg::mat &v) const { + std::hash> h; + return h(v.x); + } +}; +template +struct hash> { + std::size_t operator()(const linalg::mat &v) const { + std::hash> h; + return h(v.x) ^ (h(v.y) << M); + } +}; +template +struct hash> { + std::size_t operator()(const linalg::mat &v) const { + std::hash> h; + return h(v.x) ^ (h(v.y) << M) ^ (h(v.z) << (M * 2)); + } +}; +template +struct hash> { + std::size_t operator()(const linalg::mat &v) const { + std::hash> h; + return h(v.x) ^ (h(v.y) << M) ^ (h(v.z) << (M * 2)) ^ (h(v.w) << (M * 3)); + } +}; +} // namespace std + +// Definitions of functions too long to be defined inline +template +constexpr linalg::mat linalg::adjugate(const mat &a) { + return {{a.y.y * a.z.z - a.z.y * a.y.z, a.z.y * a.x.z - a.x.y * a.z.z, + a.x.y * a.y.z - a.y.y * a.x.z}, + {a.y.z * a.z.x - a.z.z * a.y.x, a.z.z * a.x.x - a.x.z * a.z.x, + a.x.z * a.y.x - a.y.z * a.x.x}, + {a.y.x * a.z.y - a.z.x * a.y.y, a.z.x * a.x.y - a.x.x * a.z.y, + a.x.x * a.y.y - a.y.x * a.x.y}}; +} + +template +constexpr linalg::mat linalg::adjugate(const mat &a) { + return {{a.y.y * a.z.z * a.w.w + a.w.y * a.y.z * a.z.w + + a.z.y * a.w.z * a.y.w - a.y.y * a.w.z * a.z.w - + a.z.y * a.y.z * a.w.w - a.w.y * a.z.z * a.y.w, + a.x.y * a.w.z * a.z.w + a.z.y * a.x.z * a.w.w + + a.w.y * a.z.z * a.x.w - a.w.y * a.x.z * a.z.w - + a.z.y * a.w.z * a.x.w - a.x.y * a.z.z * a.w.w, + a.x.y * a.y.z * a.w.w + a.w.y * a.x.z * a.y.w + + a.y.y * a.w.z * a.x.w - a.x.y * a.w.z * a.y.w - + a.y.y * a.x.z * a.w.w - a.w.y * a.y.z * a.x.w, + a.x.y * a.z.z * a.y.w + a.y.y * a.x.z * a.z.w + + a.z.y * a.y.z * a.x.w - a.x.y * a.y.z * a.z.w - + a.z.y * a.x.z * a.y.w - a.y.y * a.z.z * a.x.w}, + {a.y.z * a.w.w * a.z.x + a.z.z * a.y.w * a.w.x + + a.w.z * a.z.w * a.y.x - a.y.z * a.z.w * a.w.x - + a.w.z * a.y.w * a.z.x - a.z.z * a.w.w * a.y.x, + a.x.z * a.z.w * a.w.x + a.w.z * a.x.w * a.z.x + + a.z.z * a.w.w * a.x.x - a.x.z * a.w.w * a.z.x - + a.z.z * a.x.w * a.w.x - a.w.z * a.z.w * a.x.x, + a.x.z * a.w.w * a.y.x + a.y.z * a.x.w * a.w.x + + a.w.z * a.y.w * a.x.x - a.x.z * a.y.w * a.w.x - + a.w.z * a.x.w * a.y.x - a.y.z * a.w.w * a.x.x, + a.x.z * a.y.w * a.z.x + a.z.z * a.x.w * a.y.x + + a.y.z * a.z.w * a.x.x - a.x.z * a.z.w * a.y.x - + a.y.z * a.x.w * a.z.x - a.z.z * a.y.w * a.x.x}, + {a.y.w * a.z.x * a.w.y + a.w.w * a.y.x * a.z.y + + a.z.w * a.w.x * a.y.y - a.y.w * a.w.x * a.z.y - + a.z.w * a.y.x * a.w.y - a.w.w * a.z.x * a.y.y, + a.x.w * a.w.x * a.z.y + a.z.w * a.x.x * a.w.y + + a.w.w * a.z.x * a.x.y - a.x.w * a.z.x * a.w.y - + a.w.w * a.x.x * a.z.y - a.z.w * a.w.x * a.x.y, + a.x.w * a.y.x * a.w.y + a.w.w * a.x.x * a.y.y + + a.y.w * a.w.x * a.x.y - a.x.w * a.w.x * a.y.y - + a.y.w * a.x.x * a.w.y - a.w.w * a.y.x * a.x.y, + a.x.w * a.z.x * a.y.y + a.y.w * a.x.x * a.z.y + + a.z.w * a.y.x * a.x.y - a.x.w * a.y.x * a.z.y - + a.z.w * a.x.x * a.y.y - a.y.w * a.z.x * a.x.y}, + {a.y.x * a.w.y * a.z.z + a.z.x * a.y.y * a.w.z + + a.w.x * a.z.y * a.y.z - a.y.x * a.z.y * a.w.z - + a.w.x * a.y.y * a.z.z - a.z.x * a.w.y * a.y.z, + a.x.x * a.z.y * a.w.z + a.w.x * a.x.y * a.z.z + + a.z.x * a.w.y * a.x.z - a.x.x * a.w.y * a.z.z - + a.z.x * a.x.y * a.w.z - a.w.x * a.z.y * a.x.z, + a.x.x * a.w.y * a.y.z + a.y.x * a.x.y * a.w.z + + a.w.x * a.y.y * a.x.z - a.x.x * a.y.y * a.w.z - + a.w.x * a.x.y * a.y.z - a.y.x * a.w.y * a.x.z, + a.x.x * a.y.y * a.z.z + a.z.x * a.x.y * a.y.z + + a.y.x * a.z.y * a.x.z - a.x.x * a.z.y * a.y.z - + a.y.x * a.x.y * a.z.z - a.z.x * a.y.y * a.x.z}}; +} + +template +constexpr T linalg::determinant(const mat &a) { + return a.x.x * (a.y.y * a.z.z * a.w.w + a.w.y * a.y.z * a.z.w + + a.z.y * a.w.z * a.y.w - a.y.y * a.w.z * a.z.w - + a.z.y * a.y.z * a.w.w - a.w.y * a.z.z * a.y.w) + + a.x.y * (a.y.z * a.w.w * a.z.x + a.z.z * a.y.w * a.w.x + + a.w.z * a.z.w * a.y.x - a.y.z * a.z.w * a.w.x - + a.w.z * a.y.w * a.z.x - a.z.z * a.w.w * a.y.x) + + a.x.z * (a.y.w * a.z.x * a.w.y + a.w.w * a.y.x * a.z.y + + a.z.w * a.w.x * a.y.y - a.y.w * a.w.x * a.z.y - + a.z.w * a.y.x * a.w.y - a.w.w * a.z.x * a.y.y) + + a.x.w * (a.y.x * a.w.y * a.z.z + a.z.x * a.y.y * a.w.z + + a.w.x * a.z.y * a.y.z - a.y.x * a.z.y * a.w.z - + a.w.x * a.y.y * a.z.z - a.z.x * a.w.y * a.y.z); +} + +template +linalg::vec linalg::rotation_quat(const mat &m) { + const vec q{m.x.x - m.y.y - m.z.z, m.y.y - m.x.x - m.z.z, + m.z.z - m.x.x - m.y.y, m.x.x + m.y.y + m.z.z}, + s[]{{1, m.x.y + m.y.x, m.z.x + m.x.z, m.y.z - m.z.y}, + {m.x.y + m.y.x, 1, m.y.z + m.z.y, m.z.x - m.x.z}, + {m.x.z + m.z.x, m.y.z + m.z.y, 1, m.x.y - m.y.x}, + {m.y.z - m.z.y, m.z.x - m.x.z, m.x.y - m.y.x, 1}}; + return copysign(normalize(sqrt(max(T(0), T(1) + q))), s[argmax(q)]); +} + +template +linalg::mat linalg::lookat_matrix(const vec &eye, + const vec ¢er, + const vec &view_y_dir, + fwd_axis a) { + const vec f = normalize(center - eye), z = a == pos_z ? f : -f, + x = normalize(cross(view_y_dir, z)), y = cross(z, x); + return inverse(mat{{x, 0}, {y, 0}, {z, 0}, {eye, 1}}); +} + +template +linalg::mat linalg::frustum_matrix(T x0, T x1, T y0, T y1, T n, T f, + fwd_axis a, z_range z) { + const T s = a == pos_z ? T(1) : T(-1), o = z == neg_one_to_one ? n : 0; + return {{2 * n / (x1 - x0), 0, 0, 0}, + {0, 2 * n / (y1 - y0), 0, 0}, + {-s * (x0 + x1) / (x1 - x0), -s * (y0 + y1) / (y1 - y0), + s * (f + o) / (f - n), s}, + {0, 0, -(n + o) * f / (f - n), 0}}; +} + +#endif diff --git a/include/manifold/manifold.h b/include/manifold/manifold.h index 9c4b30013..0737a9ac1 100644 --- a/include/manifold/manifold.h +++ b/include/manifold/manifold.h @@ -101,22 +101,22 @@ struct MeshGLP { bool Merge(); - glm::vec<3, Precision> GetVertPos(size_t i) const { + la::vec GetVertPos(size_t i) const { size_t offset = i * numProp; - return glm::vec<3, Precision>(vertProperties[offset], - vertProperties[offset + 1], - vertProperties[offset + 2]); + return la::vec(vertProperties[offset], + vertProperties[offset + 1], + vertProperties[offset + 2]); } - glm::vec<3, I> GetTriVerts(size_t i) const { + la::vec GetTriVerts(size_t i) const { size_t offset = 3 * i; - return glm::vec<3, I>(triVerts[offset], triVerts[offset + 1], - triVerts[offset + 2]); + return la::vec(triVerts[offset], triVerts[offset + 1], + triVerts[offset + 2]); } - glm::vec<4, Precision> GetTangent(size_t i) const { + la::vec GetTangent(size_t i) const { size_t offset = 4 * i; - return glm::vec<4, Precision>( + return la::vec( halfedgeTangent[offset], halfedgeTangent[offset + 1], halfedgeTangent[offset + 2], halfedgeTangent[offset + 3]); } @@ -256,7 +256,7 @@ class Manifold { Manifold Scale(vec3) const; Manifold Rotate(double xDegrees, double yDegrees = 0.0, double zDegrees = 0.0) const; - Manifold Transform(const mat4x3&) const; + Manifold Transform(const mat3x4&) const; Manifold Mirror(vec3) const; Manifold Warp(std::function) const; Manifold WarpBatch(std::function)>) const; diff --git a/manifold.pc.in b/manifold.pc.in index a9fc25c24..edd066175 100644 --- a/manifold.pc.in +++ b/manifold.pc.in @@ -7,6 +7,6 @@ Name: manifold@PCFILE_LIB_SUFFIX@ Description: Geometry library for topological robustness Version: @MANIFOLD_VERSION@ URL: https://github.com/elalish/manifold -Requires-private: glm tbb @TEMPLATE_OPTIONAL_CLIPPER@ +Requires-private: tbb @TEMPLATE_OPTIONAL_CLIPPER@ Libs: -L${libdir} -lmanifold@PCFILE_LIB_SUFFIX@ Cflags: -I${includedir} diff --git a/manifoldConfig.cmake.in b/manifoldConfig.cmake.in index 710f0f26f..8437f5eb8 100644 --- a/manifoldConfig.cmake.in +++ b/manifoldConfig.cmake.in @@ -11,8 +11,6 @@ if(_FIND_ROOT STREQUAL "/") set(_FIND_ROOT "") endif() -set(glm_ROOT "${_FIND_ROOT}") -find_package(glm REQUIRED) set(MANIFOLD_CROSS_SECTION "@MANIFOLD_CROSS_SECTION@") if(MANIFOLD_CROSS_SECTION) set(Clipper2_ROOT "${_FIND_ROOT}") diff --git a/samples/src/bracelet.cpp b/samples/src/bracelet.cpp index 5138b0654..e97c9d642 100644 --- a/samples/src/bracelet.cpp +++ b/samples/src/bracelet.cpp @@ -35,15 +35,15 @@ Manifold Base(double width, double radius, double decorRadius, } Polygons stretch(1); - double dPhiRad = 2 * glm::pi() / nCut; + double dPhiRad = 2 * kPi / nCut; vec2 p0(outerRadius, 0.0); vec2 p1(innerRadius, -cut); vec2 p2(innerRadius, cut); for (int i = 0; i < nCut; ++i) { - stretch[0].push_back(glm::rotate(p0, dPhiRad * i)); - stretch[0].push_back(glm::rotate(p1, dPhiRad * i)); - stretch[0].push_back(glm::rotate(p2, dPhiRad * i)); - stretch[0].push_back(glm::rotate(p0, dPhiRad * i)); + stretch[0].push_back(la::rot(dPhiRad * i, p0)); + stretch[0].push_back(la::rot(dPhiRad * i, p1)); + stretch[0].push_back(la::rot(dPhiRad * i, p2)); + stretch[0].push_back(la::rot(dPhiRad * i, p0)); } base = Manifold::Extrude(stretch, width) ^ base; @@ -75,11 +75,11 @@ namespace manifold { Manifold StretchyBracelet(double radius, double height, double width, double thickness, int nDecor, int nCut, int nDivision) { - double twistRadius = glm::pi() * radius / nDecor; + double twistRadius = kPi * radius / nDecor; double decorRadius = twistRadius * 1.5; double outerRadius = radius + (decorRadius + twistRadius) * 0.5; double innerRadius = outerRadius - height; - double cut = 0.5 * (glm::pi() * 2 * innerRadius / nCut - thickness); + double cut = 0.5 * (kPi * 2 * innerRadius / nCut - thickness); double adjThickness = 0.5 * thickness * height / cut; return Base(width, radius, decorRadius, twistRadius, nDecor, diff --git a/samples/src/condensed_matter.cpp b/samples/src/condensed_matter.cpp index 3c883f823..4ae771396 100644 --- a/samples/src/condensed_matter.cpp +++ b/samples/src/condensed_matter.cpp @@ -18,7 +18,6 @@ namespace { using namespace manifold; -using namespace glm; using manifold::vec3; constexpr double AtomicRadiusN2 = 0.65; @@ -137,8 +136,9 @@ Manifold CondensedMatter(int fn) { for (int x = -3; x <= 3; x++) for (int y = -1; y <= 2; y++) parts.push_back( - GraphiteCell(fn, {x + (y % 2 == 0 ? 0.0 : 0.5), y, - LayerSeperationC * 0.5 + LatticeCellSizeSi * 1.5}) + GraphiteCell(fn, + {x + (y % 2 == 0 ? 0.0 : 0.5), static_cast(y), + LayerSeperationC * 0.5 + LatticeCellSizeSi * 1.5}) .Translate({0, -siOffset, 0}) .Rotate(0, 0, 45)); diff --git a/samples/src/gyroid_module.cpp b/samples/src/gyroid_module.cpp index 98350fe7a..486909e16 100644 --- a/samples/src/gyroid_module.cpp +++ b/samples/src/gyroid_module.cpp @@ -20,13 +20,13 @@ using namespace manifold; struct Gyroid { double operator()(vec3 p) const { - p -= glm::pi() / 4; + p -= kPi / 4; return cos(p.x) * sin(p.y) + cos(p.y) * sin(p.z) + cos(p.z) * sin(p.x); } }; Manifold RhombicDodecahedron(double size) { - Manifold box = Manifold::Cube(size * glm::sqrt(2.0) * vec3(1, 1, 2), true); + Manifold box = Manifold::Cube(size * la::sqrt(2.0) * vec3(1, 1, 2), true); Manifold result = box.Rotate(90, 45) ^ box.Rotate(90, 45, 90); return result ^ box.Rotate(0, 0, 45); } @@ -44,7 +44,7 @@ namespace manifold { */ Manifold GyroidModule(double size, int n) { auto gyroid = [&](double level) { - const double period = glm::two_pi(); + const double period = kTwoPi; return Manifold::LevelSet(Gyroid(), {vec3(-period), vec3(period)}, period / n, level) .Scale(vec3(size / period)); @@ -52,6 +52,6 @@ Manifold GyroidModule(double size, int n) { Manifold result = (RhombicDodecahedron(size) ^ gyroid(-0.4)) - gyroid(0.4); - return result.Rotate(-45, 0, 90).Translate({0, 0, size / glm::sqrt(2.0)}); + return result.Rotate(-45, 0, 90).Translate({0, 0, size / la::sqrt(2.0)}); } } // namespace manifold diff --git a/samples/src/knot.cpp b/samples/src/knot.cpp index 09d6e4350..cda74a14a 100644 --- a/samples/src/knot.cpp +++ b/samples/src/knot.cpp @@ -56,16 +56,16 @@ Manifold TorusKnot(int p, int q, double majorRadius, double minorRadius, double psi = q * atan2(v.x, v.y); double theta = psi * p / q; vec2 xy = vec2(v); - double x1 = sqrt(glm::dot(xy, xy)); + double x1 = sqrt(la::dot(xy, xy)); double phi = atan2(x1 - 2, v.z); v = vec3(cos(phi), 0.0, sin(phi)); v *= threadRadius; double r = majorRadius + minorRadius * cos(theta); - v = glm::rotateX(v, -double(atan2(p * minorRadius, q * r))); + v = la::rotx(-double(atan2(p * minorRadius, q * r)), v); v.x += minorRadius; - v = glm::rotateY(v, theta); + v = la::roty(theta, v); v.x += majorRadius; - v = glm::rotateZ(v, psi); + v = la::rotz(psi, v); }); if (kLoops > 1) { diff --git a/samples/src/scallop.cpp b/samples/src/scallop.cpp index 7e7569a3b..d8789b815 100644 --- a/samples/src/scallop.cpp +++ b/samples/src/scallop.cpp @@ -32,19 +32,19 @@ Manifold Scallop() { scallop.numProp = 3; scallop.vertProperties = {-offset, 0, height, -offset, 0, -height}; - const double delta = glm::pi() / wiggles; + const double delta = kPi / wiggles; for (int i = 0; i < 2 * wiggles; ++i) { double theta = (i - wiggles) * delta; - double amp = 0.5 * height * glm::max(glm::cos(0.8 * theta), 0.0); + double amp = 0.5 * height * la::max(la::cos(0.8 * theta), 0.0); scallop.vertProperties.insert( scallop.vertProperties.end(), - {radius * glm::cos(theta), radius * glm::sin(theta), + {radius * la::cos(theta), radius * la::sin(theta), amp * (i % 2 == 0 ? 1 : -1)}); int j = i + 1; if (j == 2 * wiggles) j = 0; - double smoothness = 1 - sharpness * glm::cos((theta + delta / 2) / 2); + double smoothness = 1 - sharpness * la::cos((theta + delta / 2) / 2); size_t halfedge = scallop.triVerts.size() + 1; sharpenedEdges.push_back({halfedge, smoothness}); scallop.triVerts.insert( diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 70b2293a2..63ba4ba66 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -60,7 +60,6 @@ if(MANIFOLD_EXPORT) endif() # Dependency libraries -set(MANIFOLD_LIBS glm::glm) if(MANIFOLD_CROSS_SECTION) if(TARGET Clipper2::Clipper2) list(APPEND MANIFOLD_LIBS Clipper2::Clipper2) diff --git a/src/boolean3.cpp b/src/boolean3.cpp index d2bd3825a..4b37eee8f 100644 --- a/src/boolean3.cpp +++ b/src/boolean3.cpp @@ -123,9 +123,9 @@ inline std::pair Shadow01( yz01 = Interpolate(vertPosQ[q1s], vertPosQ[q1e], vertPosP[p0].x); if (reverse) { vec3 diff = vertPosQ[q1s] - vertPosP[p0]; - const double start2 = glm::dot(diff, diff); + const double start2 = la::dot(diff, diff); diff = vertPosQ[q1e] - vertPosP[p0]; - const double end2 = glm::dot(diff, diff); + const double end2 = la::dot(diff, diff); const double dir = start2 < end2 ? normalP[q1s].y : normalP[q1e].y; if (!Shadows(yz01[0], vertPosP[p0].y, expandP * dir)) s01 = 0; } else { @@ -211,7 +211,7 @@ struct Kernel11 { if (k < 2 && (k == 0 || (s01 != 0) != shadows)) { shadows = s01 != 0; pRL[k] = vertPosP[p0[i]]; - qRL[k] = vec3(pRL[k].x, yz01); + qRL[k] = vec3(pRL[k].x, yz01.x, yz01.y); ++k; } } @@ -229,7 +229,7 @@ struct Kernel11 { if (k < 2 && (k == 0 || (s10 != 0) != shadows)) { shadows = s10 != 0; qRL[k] = vertPosQ[q0[i]]; - pRL[k] = vec3(qRL[k].x, yz10); + pRL[k] = vec3(qRL[k].x, yz10.x, yz10.y); ++k; } } @@ -244,9 +244,9 @@ struct Kernel11 { const int p1s = halfedgeP[p1].startVert; const int p1e = halfedgeP[p1].endVert; vec3 diff = vertPosP[p1s] - vec3(xyzz11); - const double start2 = glm::dot(diff, diff); + const double start2 = la::dot(diff, diff); diff = vertPosP[p1e] - vec3(xyzz11); - const double end2 = glm::dot(diff, diff); + const double end2 = la::dot(diff, diff); const double dir = start2 < end2 ? normalP[p1s].z : normalP[p1e].z; if (!Shadows(xyzz11.z, xyzz11.w, expandP * dir)) s11 = 0; @@ -307,7 +307,7 @@ struct Kernel02 { if (!forward) { const int qVert = halfedgeQ[q1F].startVert; const vec3 diff = posP - vertPosQ[qVert]; - const double metric = glm::dot(diff, diff); + const double metric = la::dot(diff, diff); if (metric < minMetric) { minMetric = metric; closestVert = qVert; diff --git a/src/boolean_result.cpp b/src/boolean_result.cpp index d90578248..4e68f1b24 100644 --- a/src/boolean_result.cpp +++ b/src/boolean_result.cpp @@ -320,12 +320,12 @@ void AppendPartialEdges(Manifold::Impl &outR, Vec &wholeHalfedgeP, const vec3 edgeVec = vertPosP[vEnd] - vertPosP[vStart]; // Fill in the edge positions of the old points. for (EdgePos &edge : edgePosP) { - edge.edgePos = glm::dot(outR.vertPos_[edge.vert], edgeVec); + edge.edgePos = la::dot(outR.vertPos_[edge.vert], edgeVec); } int inclusion = i03[vStart]; EdgePos edgePos = {vP2R[vStart], - glm::dot(outR.vertPos_[vP2R[vStart]], edgeVec), + la::dot(outR.vertPos_[vP2R[vStart]], edgeVec), inclusion > 0}; for (int j = 0; j < std::abs(inclusion); ++j) { edgePosP.push_back(edgePos); @@ -333,7 +333,7 @@ void AppendPartialEdges(Manifold::Impl &outR, Vec &wholeHalfedgeP, } inclusion = i03[vEnd]; - edgePos = {vP2R[vEnd], glm::dot(outR.vertPos_[vP2R[vEnd]], edgeVec), + edgePos = {vP2R[vEnd], la::dot(outR.vertPos_[vP2R[vEnd]], edgeVec), inclusion < 0}; for (int j = 0; j < std::abs(inclusion); ++j) { edgePosP.push_back(edgePos); @@ -645,7 +645,7 @@ void CreateProperties(Manifold::Impl &outR, const Manifold::Impl &inP, vec3 oldProps; for (const int j : {0, 1, 2}) oldProps[j] = properties[oldNumProp * triProp[j] + p]; - outR.meshRelation_.properties.push_back(glm::dot(uvw, oldProps)); + outR.meshRelation_.properties.push_back(la::dot(uvw, oldProps)); } else { outR.meshRelation_.properties.push_back(0); } diff --git a/src/collider.h b/src/collider.h index 5cc570408..6c303fb5c 100644 --- a/src/collider.h +++ b/src/collider.h @@ -253,7 +253,7 @@ struct BuildInternalBoxes { }; struct TransformBox { - const mat4x3 transform; + const mat3x4 transform; void operator()(Box& box) { box = box.Transform(transform); } }; @@ -288,7 +288,7 @@ class Collider { UpdateBoxes(leafBB); } - bool Transform(mat4x3 transform) { + bool Transform(mat3x4 transform) { ZoneScoped; bool axisAligned = true; for (int row : {0, 1, 2}) { @@ -360,7 +360,7 @@ class Collider { static uint32_t MortonCode(vec3 position, Box bBox) { using collider_internal::SpreadBits3; vec3 xyz = (position - bBox.min) / (bBox.max - bBox.min); - xyz = glm::min(vec3(1023.0), glm::max(vec3(0.0), 1024.0 * xyz)); + xyz = la::min(vec3(1023.0), la::max(vec3(0.0), 1024.0 * xyz)); uint32_t x = SpreadBits3(static_cast(xyz.x)); uint32_t y = SpreadBits3(static_cast(xyz.y)); uint32_t z = SpreadBits3(static_cast(xyz.z)); diff --git a/src/constructors.cpp b/src/constructors.cpp index 594066c22..34b89eddd 100644 --- a/src/constructors.cpp +++ b/src/constructors.cpp @@ -113,10 +113,11 @@ Manifold Manifold::Tetrahedron() { * @param center Set to true to shift the center to the origin. */ Manifold Manifold::Cube(vec3 size, bool center) { - if (size.x < 0.0 || size.y < 0.0 || size.z < 0.0 || glm::length(size) == 0.) { + if (size.x < 0.0 || size.y < 0.0 || size.z < 0.0 || la::length(size) == 0.) { return Invalid(); } - mat4x3 m(glm::translate(center ? (-size / 2.0) : vec3(0)) * glm::scale(size)); + mat3x4 m({{size.x, 0.0, 0.0}, {0.0, size.y, 0.0}, {0.0, 0.0, size.z}}, + center ? (-size / 2.0) : vec3(0.0)); return Manifold(std::make_shared(Manifold::Impl::Shape::Cube, m)); } @@ -176,8 +177,8 @@ Manifold Manifold::Sphere(double radius, int circularSegments) { [n](vec3 edge, vec4 tangentStart, vec4 tangentEnd) { return n - 1; }); for_each_n(autoPolicy(pImpl_->NumVert(), 1e5), pImpl_->vertPos_.begin(), pImpl_->NumVert(), [radius](vec3& v) { - v = glm::cos(glm::half_pi() * (1.0 - v)); - v = radius * glm::normalize(v); + v = la::cos(kHalfPi * (1.0 - v)); + v = radius * la::normalize(v); if (std::isnan(v.x)) v = vec3(0.0); }); pImpl_->Finish(); @@ -236,9 +237,9 @@ Manifold Manifold::Extrude(const Polygons& crossSection, double height, for (int i = 1; i < nDivisions + 1; ++i) { double alpha = i / double(nDivisions); double phi = alpha * twistDegrees; - vec2 scale = glm::mix(vec2(1.0), scaleTop, alpha); - mat2 rotation(cosd(phi), sind(phi), -sind(phi), cosd(phi)); - mat2 transform = mat2(scale.x, 0.0, 0.0, scale.y) * rotation; + vec2 scale = la::lerp(vec2(1.0), scaleTop, alpha); + mat2 rotation({cosd(phi), sind(phi)}, {-sind(phi), cosd(phi)}); + mat2 transform = mat2({scale.x, 0.0}, {0.0, scale.y}) * rotation; size_t j = 0; size_t idx = 0; for (const auto& poly : crossSection) { @@ -247,14 +248,16 @@ Manifold Manifold::Extrude(const Polygons& crossSection, double height, size_t thisVert = vert + offset; size_t lastVert = (vert == 0 ? poly.size() : vert) - 1 + offset; if (i == nDivisions && isCone) { - triVerts.push_back({nCrossSection * i + j, lastVert - nCrossSection, - thisVert - nCrossSection}); + triVerts.push_back(ivec3(nCrossSection * i + j, + lastVert - nCrossSection, + thisVert - nCrossSection)); } else { vec2 pos = transform * poly[vert]; vertPos.push_back({pos.x, pos.y, height * alpha}); - triVerts.push_back({thisVert, lastVert, thisVert - nCrossSection}); triVerts.push_back( - {lastVert, lastVert - nCrossSection, thisVert - nCrossSection}); + ivec3(thisVert, lastVert, thisVert - nCrossSection)); + triVerts.push_back(ivec3(lastVert, lastVert - nCrossSection, + thisVert - nCrossSection)); } } ++j; @@ -383,18 +386,18 @@ Manifold Manifold::Revolve(const Polygons& crossSection, int circularSegments, if (isFullRevolution || slice > 0) { const int lastSlice = (slice == 0 ? nDivisions : slice) - 1; if (currPolyVertex.x > 0.0) { - triVerts.push_back( - {startPosIndex + slice, startPosIndex + lastSlice, - // "Reuse" vertex of first slice if it lies on the revolve axis - (prevPolyVertex.x == 0.0 ? prevStartPosIndex - : prevStartPosIndex + lastSlice)}); + triVerts.push_back(ivec3( + startPosIndex + slice, startPosIndex + lastSlice, + // "Reuse" vertex of first slice if it lies on the revolve axis + (prevPolyVertex.x == 0.0 ? prevStartPosIndex + : prevStartPosIndex + lastSlice))); } if (prevPolyVertex.x > 0.0) { triVerts.push_back( - {prevStartPosIndex + lastSlice, prevStartPosIndex + slice, - (currPolyVertex.x == 0.0 ? startPosIndex - : startPosIndex + slice)}); + ivec3(prevStartPosIndex + lastSlice, prevStartPosIndex + slice, + (currPolyVertex.x == 0.0 ? startPosIndex + : startPosIndex + slice))); } } } diff --git a/src/cross_section/cross_section.cpp b/src/cross_section/cross_section.cpp index 6e78331b1..f46013433 100644 --- a/src/cross_section/cross_section.cpp +++ b/src/cross_section/cross_section.cpp @@ -95,8 +95,8 @@ C2::PathD pathd_of_contour(const SimplePolygon& ctr) { return p; } -C2::PathsD transform(const C2::PathsD ps, const mat3x2 m) { - const bool invert = glm::determinant(mat2(m)) < 0; +C2::PathsD transform(const C2::PathsD ps, const mat2x3 m) { + const bool invert = la::determinant(mat2(m)) < 0; auto transformed = C2::PathsD(); transformed.reserve(ps.size()); for (auto path : ps) { @@ -292,11 +292,11 @@ CrossSection::CrossSection(const Rect& rect) { // All access to paths_ should be done through the GetPaths() method, which // applies the accumulated transform_ std::shared_ptr CrossSection::GetPaths() const { - if (transform_ == mat3x2(1.0)) { + if (transform_ == Identity2x3()) { return paths_; } paths_ = shared_paths(::transform(paths_->paths_, transform_)); - transform_ = mat3x2(1.0); + transform_ = Identity2x3(); return paths_; } @@ -309,7 +309,7 @@ std::shared_ptr CrossSection::GetPaths() const { * @param center Set to true to shift the center to the origin. */ CrossSection CrossSection::Square(const vec2 size, bool center) { - if (size.x < 0.0 || size.y < 0.0 || glm::length(size) == 0.0) { + if (size.x < 0.0 || size.y < 0.0 || la::length(size) == 0.0) { return CrossSection(); } @@ -483,9 +483,9 @@ std::vector CrossSection::Decompose() const { * @param v The vector to add to every vertex. */ CrossSection CrossSection::Translate(const vec2 v) const { - mat3x2 m(1.0, 0.0, // - 0.0, 1.0, // - v.x, v.y); + mat2x3 m({1.0, 0.0}, // + {0.0, 1.0}, // + {v.x, v.y}); return Transform(m); } @@ -498,9 +498,9 @@ CrossSection CrossSection::Translate(const vec2 v) const { CrossSection CrossSection::Rotate(double degrees) const { auto s = sind(degrees); auto c = cosd(degrees); - mat3x2 m(c, s, // - -s, c, // - 0.0, 0.0); + mat2x3 m({c, s}, // + {-s, c}, // + {0.0, 0.0}); return Transform(m); } @@ -511,9 +511,9 @@ CrossSection CrossSection::Rotate(double degrees) const { * @param v The vector to multiply every vertex by per component. */ CrossSection CrossSection::Scale(const vec2 scale) const { - mat3x2 m(scale.x, 0.0, // - 0.0, scale.y, // - 0.0, 0.0); + mat2x3 m({scale.x, 0.0}, // + {0.0, scale.y}, // + {0.0, 0.0}); return Transform(m); } @@ -526,11 +526,11 @@ CrossSection CrossSection::Scale(const vec2 scale) const { * @param ax the axis to be mirrored over */ CrossSection CrossSection::Mirror(const vec2 ax) const { - if (glm::length(ax) == 0.) { + if (la::length(ax) == 0.) { return CrossSection(); } - auto n = glm::normalize(glm::abs(ax)); - auto m = mat3x2(mat2(1.0) - 2.0 * glm::outerProduct(n, n)); + auto n = la::normalize(la::abs(ax)); + auto m = mat2x3(mat2(la::identity) - 2.0 * la::outerprod(n, n), vec2(0.0)); return Transform(m); } @@ -541,9 +541,9 @@ CrossSection CrossSection::Mirror(const vec2 ax) const { * * @param m The affine transform matrix to apply to all the vertices. */ -CrossSection CrossSection::Transform(const mat3x2& m) const { +CrossSection CrossSection::Transform(const mat2x3& m) const { auto transformed = CrossSection(); - transformed.transform_ = m * mat3(transform_); + transformed.transform_ = m * Mat3(transform_); transformed.paths_ = paths_; return transformed; } diff --git a/src/csg_tree.cpp b/src/csg_tree.cpp index 2b99f8b0f..3e6cbe777 100644 --- a/src/csg_tree.cpp +++ b/src/csg_tree.cpp @@ -32,7 +32,7 @@ constexpr int kParallelThreshold = 4096; namespace { using namespace manifold; struct Transform4x3 { - mat4x3 transform; + mat3x4 transform; vec3 operator()(vec3 position) const { return transform * vec4(position, 1.0); @@ -111,29 +111,29 @@ std::shared_ptr CsgNode::Boolean( } std::shared_ptr CsgNode::Translate(const vec3 &t) const { - mat4x3 transform(1.0); + mat3x4 transform = Identity3x4(); transform[3] += t; return Transform(transform); } std::shared_ptr CsgNode::Scale(const vec3 &v) const { - mat4x3 transform(1.0); - for (int i : {0, 1, 2}) transform[i] *= v; + mat3x4 transform; + for (int i : {0, 1, 2}) transform[i][i] = v[i]; return Transform(transform); } std::shared_ptr CsgNode::Rotate(double xDegrees, double yDegrees, double zDegrees) const { - mat3 rX(1.0, 0.0, 0.0, // - 0.0, cosd(xDegrees), sind(xDegrees), // - 0.0, -sind(xDegrees), cosd(xDegrees)); - mat3 rY(cosd(yDegrees), 0.0, -sind(yDegrees), // - 0.0, 1.0, 0.0, // - sind(yDegrees), 0.0, cosd(yDegrees)); - mat3 rZ(cosd(zDegrees), sind(zDegrees), 0.0, // - -sind(zDegrees), cosd(zDegrees), 0.0, // - 0.0, 0.0, 1.0); - mat4x3 transform(rZ * rY * rX); + mat3 rX({1.0, 0.0, 0.0}, // + {0.0, cosd(xDegrees), sind(xDegrees)}, // + {0.0, -sind(xDegrees), cosd(xDegrees)}); + mat3 rY({cosd(yDegrees), 0.0, -sind(yDegrees)}, // + {0.0, 1.0, 0.0}, // + {sind(yDegrees), 0.0, cosd(yDegrees)}); + mat3 rZ({cosd(zDegrees), sind(zDegrees), 0.0}, // + {-sind(zDegrees), cosd(zDegrees), 0.0}, // + {0.0, 0.0, 1.0}); + mat3x4 transform(rZ * rY * rX, vec3()); return Transform(transform); } @@ -143,25 +143,25 @@ CsgLeafNode::CsgLeafNode(std::shared_ptr pImpl_) : pImpl_(pImpl_) {} CsgLeafNode::CsgLeafNode(std::shared_ptr pImpl_, - mat4x3 transform_) + mat3x4 transform_) : pImpl_(pImpl_), transform_(transform_) {} std::shared_ptr CsgLeafNode::GetImpl() const { - if (transform_ == mat4x3(1.0)) return pImpl_; + if (transform_ == Identity3x4()) return pImpl_; pImpl_ = std::make_shared(pImpl_->Transform(transform_)); - transform_ = mat4x3(1.0); + transform_ = Identity3x4(); return pImpl_; } -mat4x3 CsgLeafNode::GetTransform() const { return transform_; } +mat3x4 CsgLeafNode::GetTransform() const { return transform_; } std::shared_ptr CsgLeafNode::ToLeafNode() const { return std::make_shared(*this); } -std::shared_ptr CsgLeafNode::Transform(const mat4x3 &m) const { - return std::make_shared(pImpl_, m * mat4(transform_)); +std::shared_ptr CsgLeafNode::Transform(const mat3x4 &m) const { + return std::make_shared(pImpl_, m * Mat4(transform_)); } CsgNodeType CsgLeafNode::GetNodeType() const { return CsgNodeType::Leaf; } @@ -270,7 +270,7 @@ Manifold::Impl CsgLeafNode::Compose( } } - if (node->transform_ == mat4x3(1.0)) { + if (node->transform_ == Identity3x4()) { copy(node->pImpl_->vertPos_.begin(), node->pImpl_->vertPos_.end(), combined.vertPos_.begin() + vertIndices[i]); copy(node->pImpl_->faceNormal_.begin(), @@ -282,7 +282,7 @@ Manifold::Impl CsgLeafNode::Compose( auto vertPosBegin = TransformIterator( node->pImpl_->vertPos_.begin(), Transform4x3({node->transform_})); mat3 normalTransform = - glm::inverse(glm::transpose(mat3(node->transform_))); + la::inverse(la::transpose(mat3(node->transform_))); auto faceNormalBegin = TransformIterator(node->pImpl_->faceNormal_.begin(), TransformNormals({normalTransform})); @@ -291,7 +291,7 @@ Manifold::Impl CsgLeafNode::Compose( copy_n(faceNormalBegin, node->pImpl_->faceNormal_.size(), combined.faceNormal_.begin() + triIndices[i]); - const bool invert = glm::determinant(mat3(node->transform_)) < 0; + const bool invert = la::determinant(mat3(node->transform_)) < 0; for_each_n(policy, countAt(0), node->pImpl_->halfedgeTangent_.size(), TransformTangents{combined.halfedgeTangent_, edgeIndices[i], mat3(node->transform_), @@ -351,14 +351,13 @@ std::shared_ptr CsgOpNode::Boolean( auto isReused = [](const auto &node) { return node->impl_.UseCount() > 1; }; - auto copyChildren = [&](const auto &list, const mat4x3 &transform) { + auto copyChildren = [&](const auto &list, const mat3x4 &transform) { for (const auto &child : list) { children.push_back(child->Transform(transform)); } }; auto self = std::dynamic_pointer_cast(shared_from_this()); - assert(self); if (IsOp(op) && !isReused(self)) { auto impl = impl_.GetGuard(); copyChildren(impl->children_, transform_); @@ -389,10 +388,10 @@ std::shared_ptr CsgOpNode::Boolean( return std::make_shared(children, op); } -std::shared_ptr CsgOpNode::Transform(const mat4x3 &m) const { +std::shared_ptr CsgOpNode::Transform(const mat3x4 &m) const { auto node = std::make_shared(); node->impl_ = impl_; - node->transform_ = m * mat4(transform_); + node->transform_ = m * Mat4(transform_); node->op_ = op_; return node; } @@ -641,6 +640,6 @@ bool CsgOpNode::IsOp(OpType op) { } } -mat4x3 CsgOpNode::GetTransform() const { return transform_; } +mat3x4 CsgOpNode::GetTransform() const { return transform_; } } // namespace manifold diff --git a/src/csg_tree.h b/src/csg_tree.h index 5e318dc6c..37c7025c0 100644 --- a/src/csg_tree.h +++ b/src/csg_tree.h @@ -25,9 +25,9 @@ class CsgLeafNode; class CsgNode : public std::enable_shared_from_this { public: virtual std::shared_ptr ToLeafNode() const = 0; - virtual std::shared_ptr Transform(const mat4x3 &m) const = 0; + virtual std::shared_ptr Transform(const mat3x4 &m) const = 0; virtual CsgNodeType GetNodeType() const = 0; - virtual mat4x3 GetTransform() const = 0; + virtual mat3x4 GetTransform() const = 0; virtual std::shared_ptr Boolean( const std::shared_ptr &second, OpType op); @@ -42,24 +42,24 @@ class CsgLeafNode final : public CsgNode { public: CsgLeafNode(); CsgLeafNode(std::shared_ptr pImpl_); - CsgLeafNode(std::shared_ptr pImpl_, mat4x3 transform_); + CsgLeafNode(std::shared_ptr pImpl_, mat3x4 transform_); std::shared_ptr GetImpl() const; std::shared_ptr ToLeafNode() const override; - std::shared_ptr Transform(const mat4x3 &m) const override; + std::shared_ptr Transform(const mat3x4 &m) const override; CsgNodeType GetNodeType() const override; - mat4x3 GetTransform() const override; + mat3x4 GetTransform() const override; static Manifold::Impl Compose( const std::vector> &nodes); private: mutable std::shared_ptr pImpl_; - mutable mat4x3 transform_ = mat4x3(1.0); + mutable mat3x4 transform_ = Identity3x4(); }; class CsgOpNode final : public CsgNode { @@ -73,13 +73,13 @@ class CsgOpNode final : public CsgNode { std::shared_ptr Boolean(const std::shared_ptr &second, OpType op) override; - std::shared_ptr Transform(const mat4x3 &m) const override; + std::shared_ptr Transform(const mat3x4 &m) const override; std::shared_ptr ToLeafNode() const override; CsgNodeType GetNodeType() const override { return op_; } - mat4x3 GetTransform() const override; + mat3x4 GetTransform() const override; private: struct Impl { @@ -88,7 +88,7 @@ class CsgOpNode final : public CsgNode { }; mutable ConcurrentSharedPtr impl_ = ConcurrentSharedPtr(Impl{}); CsgNodeType op_; - mat4x3 transform_ = mat4x3(1.0); + mat3x4 transform_ = Identity3x4(); // the following fields are for lazy evaluation, so they are mutable mutable std::shared_ptr cache_ = nullptr; diff --git a/src/edge_op.cpp b/src/edge_op.cpp index 759a81031..9398deade 100644 --- a/src/edge_op.cpp +++ b/src/edge_op.cpp @@ -29,7 +29,7 @@ ivec3 TriOf(int edge) { bool Is01Longest(vec2 v0, vec2 v1, vec2 v2) { const vec2 e[3] = {v1 - v0, v2 - v1, v0 - v2}; double l[3]; - for (int i : {0, 1, 2}) l[i] = glm::dot(e[i], e[i]); + for (int i : {0, 1, 2}) l[i] = la::dot(e[i], e[i]); return l[0] > l[1] && l[0] > l[2]; } @@ -54,7 +54,7 @@ struct ShortEdge { // Flag short edges const vec3 delta = vertPos[halfedge[edge].endVert] - vertPos[halfedge[edge].startVert]; - return glm::dot(delta, delta) < precision * precision; + return la::dot(delta, delta) < precision * precision; } }; @@ -90,7 +90,7 @@ struct SwappableEdge { int tri = edge / 3; ivec3 triEdge = TriOf(edge); - mat3x2 projection = GetAxisAlignedProjection(triNormal[tri]); + mat2x3 projection = GetAxisAlignedProjection(triNormal[tri]); vec2 v[3]; for (int i : {0, 1, 2}) v[i] = projection * vertPos[halfedge[triEdge[i]].startVert]; @@ -470,7 +470,7 @@ void Manifold::Impl::CollapseEdge(const int edge, std::vector& edges) { const vec3 pNew = vertPos_[endVert]; const vec3 pOld = vertPos_[toRemove.startVert]; const vec3 delta = pNew - pOld; - const bool shortEdge = glm::dot(delta, delta) < precision_ * precision_; + const bool shortEdge = la::dot(delta, delta) < precision_ * precision_; // Orbit endVert int current = halfedge_[tri0edge[1]].pairedHalfedge; @@ -491,7 +491,7 @@ void Manifold::Impl::CollapseEdge(const int edge, std::vector& edges) { vec3 pNext = vertPos_[halfedge_[current].endVert]; const int tri = current / 3; const TriRef ref = triRef[tri]; - const mat3x2 projection = GetAxisAlignedProjection(faceNormal_[tri]); + const mat2x3 projection = GetAxisAlignedProjection(faceNormal_[tri]); // Don't collapse if the edge is not redundant (this may have changed due // to the collapse of neighbors). if (!ref.SameFace(refCheck)) { @@ -577,7 +577,7 @@ void Manifold::Impl::RecursiveEdgeSwap(const int edge, int& tag, const ivec3 perm0 = TriOf(edge % 3); const ivec3 perm1 = TriOf(pair % 3); - mat3x2 projection = GetAxisAlignedProjection(faceNormal_[edge / 3]); + mat2x3 projection = GetAxisAlignedProjection(faceNormal_[edge / 3]); vec2 v[4]; for (int i : {0, 1, 2}) v[i] = projection * vertPos_[halfedge_[tri0edge[i]].startVert]; @@ -607,8 +607,8 @@ void Manifold::Impl::RecursiveEdgeSwap(const int edge, int& tag, const int tri1 = tri1edge[0] / 3; faceNormal_[tri0] = faceNormal_[tri1]; triRef[tri0] = triRef[tri1]; - const double l01 = glm::length(v[1] - v[0]); - const double l02 = glm::length(v[2] - v[0]); + const double l01 = la::length(v[1] - v[0]); + const double l02 = la::length(v[2] - v[0]); const double a = std::max(0.0, std::min(1.0, l02 / l01)); // Update properties if applicable if (meshRelation_.properties.size() > 0) { @@ -649,7 +649,7 @@ void Manifold::Impl::RecursiveEdgeSwap(const int edge, int& tag, // Two facing, long-edge degenerates can swap. SwapEdge(); const vec2 e23 = v[3] - v[2]; - if (glm::dot(e23, e23) < precision_ * precision_) { + if (la::dot(e23, e23) < precision_ * precision_) { tag++; CollapseEdge(tri0edge[2], edges); edges.resize(0); diff --git a/src/face_op.cpp b/src/face_op.cpp index b626a0117..eb65d278e 100644 --- a/src/face_op.cpp +++ b/src/face_op.cpp @@ -75,7 +75,7 @@ void Manifold::Impl::Face2Tri(const Vec& faceEdge, halfedge_[firstEdge + 1].startVert, halfedge_[firstEdge + 2].startVert, halfedge_[firstEdge + 3].startVert}; - const mat3x2 projection = GetAxisAlignedProjection(normal); + const mat2x3 projection = GetAxisAlignedProjection(normal); auto triCCW = [&projection, this](const ivec3 tri) { return CCW(projection * this->vertPos_[tri[0]], projection * this->vertPos_[tri[1]], @@ -94,8 +94,8 @@ void Manifold::Impl::Face2Tri(const Vec& faceEdge, tri1[1] = halfedge_[firstEdge + i].startVert; } } - DEBUG_ASSERT(glm::all(glm::greaterThanEqual(tri0, ivec3(0))) && - glm::all(glm::greaterThanEqual(tri1, ivec3(0))), + DEBUG_ASSERT(la::all(la::gequal(tri0, ivec3(0))) && + la::all(la::gequal(tri1, ivec3(0))), topologyErr, "non-manifold quad!"); bool firstValid = triCCW(tri0) && triCCW(tri1); tri0[2] = tri1[1]; @@ -108,8 +108,8 @@ void Manifold::Impl::Face2Tri(const Vec& faceEdge, } else if (firstValid) { vec3 firstCross = vertPos_[tri0[0]] - vertPos_[tri1[0]]; vec3 secondCross = vertPos_[tri0[1]] - vertPos_[tri1[1]]; - if (glm::dot(firstCross, firstCross) < - glm::dot(secondCross, secondCross)) { + if (la::dot(firstCross, firstCross) < + la::dot(secondCross, secondCross)) { tri0[2] = tri1[0]; tri1[2] = tri0[0]; } @@ -126,7 +126,7 @@ void Manifold::Impl::Face2Tri(const Vec& faceEdge, }; auto generalTriangulation = [&](int face) { const vec3 normal = faceNormal_[face]; - const mat3x2 projection = GetAxisAlignedProjection(normal); + const mat2x3 projection = GetAxisAlignedProjection(normal); const PolygonsIdx polys = Face2Polygons(halfedge_.cbegin() + faceEdge[face], halfedge_.cbegin() + faceEdge[face + 1], projection); @@ -199,7 +199,7 @@ void Manifold::Impl::Face2Tri(const Vec& faceEdge, */ PolygonsIdx Manifold::Impl::Face2Polygons(VecView::IterC start, VecView::IterC end, - mat3x2 projection) const { + mat2x3 projection) const { std::multimap vert_edge; for (auto edge = start; edge != end; ++edge) { vert_edge.emplace( @@ -275,7 +275,7 @@ Polygons Manifold::Impl::Slice(double height) const { const vec3 below = vertPos_[up.startVert]; const vec3 above = vertPos_[up.endVert]; const double a = (height - below.z) / (above.z - below.z); - poly.push_back(vec2(glm::mix(below, above, a))); + poly.push_back(vec2(la::lerp(below, above, a))); const int pair = up.pairedHalfedge; tri = pair / 3; @@ -289,7 +289,7 @@ Polygons Manifold::Impl::Slice(double height) const { } Polygons Manifold::Impl::Project() const { - const mat3x2 projection = GetAxisAlignedProjection({0, 0, 1}); + const mat2x3 projection = GetAxisAlignedProjection({0, 0, 1}); Vec cusps(NumEdge()); cusps.resize( copy_if( diff --git a/src/impl.cpp b/src/impl.cpp index abcf6f375..9857d92bd 100644 --- a/src/impl.cpp +++ b/src/impl.cpp @@ -40,7 +40,7 @@ void AtomicAddVec3(vec3& target, const vec3& add) { } struct Transform4x3 { - const mat4x3 transform; + const mat3x4 transform; vec3 operator()(vec3 position) { return transform * vec4(position, 1.0); } }; @@ -62,21 +62,21 @@ struct AssignNormals { vec3 edge[3]; for (int i : {0, 1, 2}) { const int j = (i + 1) % 3; - edge[i] = glm::normalize(vertPos[triVerts[j]] - vertPos[triVerts[i]]); + edge[i] = la::normalize(vertPos[triVerts[j]] - vertPos[triVerts[i]]); } if (calculateTriNormal) { - triNormal = glm::normalize(glm::cross(edge[0], edge[1])); + triNormal = la::normalize(la::cross(edge[0], edge[1])); if (std::isnan(triNormal.x)) triNormal = vec3(0, 0, 1); } // corner angles vec3 phi; - double dot = -glm::dot(edge[2], edge[0]); - phi[0] = dot >= 1 ? 0 : (dot <= -1 ? glm::pi() : std::acos(dot)); - dot = -glm::dot(edge[0], edge[1]); - phi[1] = dot >= 1 ? 0 : (dot <= -1 ? glm::pi() : std::acos(dot)); - phi[2] = glm::pi() - phi[0] - phi[1]; + double dot = -la::dot(edge[2], edge[0]); + phi[0] = dot >= 1 ? 0 : (dot <= -1 ? kPi : std::acos(dot)); + dot = -la::dot(edge[0], edge[1]); + phi[1] = dot >= 1 ? 0 : (dot <= -1 ? kPi : std::acos(dot)); + phi[2] = kPi - phi[0] - phi[1]; // assign weighted sum for (int i : {0, 1, 2}) { @@ -143,18 +143,18 @@ struct CoplanarEdge { const vec3 pairVec = vertPos[halfedge[3 * pairFace + pairNum].startVert] - base; - const double length = std::max(glm::length(jointVec), glm::length(edgeVec)); + const double length = std::max(la::length(jointVec), la::length(edgeVec)); const double lengthPair = - std::max(glm::length(jointVec), glm::length(pairVec)); - vec3 normal = glm::cross(jointVec, edgeVec); - const double area = glm::length(normal); - const double areaPair = glm::length(glm::cross(pairVec, jointVec)); + std::max(la::length(jointVec), la::length(pairVec)); + vec3 normal = la::cross(jointVec, edgeVec); + const double area = la::length(normal); + const double areaPair = la::length(la::cross(pairVec, jointVec)); triArea[edgeFace] = area; triArea[pairFace] = areaPair; // Don't link degenerate triangles if (area < length * precision || areaPair < lengthPair * precision) return; - const double volume = std::abs(glm::dot(normal, pairVec)); + const double volume = std::abs(la::dot(normal, pairVec)); // Only operate on coplanar triangles if (volume > std::max(area, areaPair) * precision) return; @@ -175,10 +175,10 @@ struct CoplanarEdge { const vec3 iEdgeVec = edgeVec + normal * scale * (edgeProp - baseProp); const vec3 iPairVec = pairVec + normal * scale * (pairProp - baseProp); - vec3 cross = glm::cross(iJointVec, iEdgeVec); + vec3 cross = la::cross(iJointVec, iEdgeVec); const double areaP = std::max( - glm::length(cross), glm::length(glm::cross(iPairVec, iJointVec))); - const double volumeP = std::abs(glm::dot(cross, iPairVec)); + la::length(cross), la::length(la::cross(iPairVec, iJointVec))); + const double volumeP = std::abs(la::dot(cross, iPairVec)); // Only operate on consistent triangles if (volumeP > areaP * precision) return; } @@ -201,16 +201,16 @@ struct CheckCoplanarity { if (referenceTri < 0 || referenceTri == tri) return; const vec3 origin = vertPos[halfedge[3 * referenceTri].startVert]; - const vec3 normal = glm::normalize( - glm::cross(vertPos[halfedge[3 * referenceTri + 1].startVert] - origin, - vertPos[halfedge[3 * referenceTri + 2].startVert] - origin)); + const vec3 normal = la::normalize( + la::cross(vertPos[halfedge[3 * referenceTri + 1].startVert] - origin, + vertPos[halfedge[3 * referenceTri + 2].startVert] - origin)); for (const int i : {0, 1, 2}) { const vec3 vert = vertPos[halfedge[3 * tri + i].startVert]; // If any component vertex is not coplanar with the component's reference // triangle, unmark the entire component so that none of its triangles are // marked coplanar. - if (std::abs(glm::dot(normal, vert - origin)) > precision) { + if (std::abs(la::dot(normal, vert - origin)) > precision) { reinterpret_cast*>(&comp2tri[component]) ->store(-1, std::memory_order_relaxed); break; @@ -255,7 +255,7 @@ uint32_t Manifold::Impl::ReserveIDs(uint32_t n) { * Create either a unit tetrahedron, cube or octahedron. The cube is in the * first octant, while the others are symmetric about the origin. */ -Manifold::Impl::Impl(Shape shape, const mat4x3 m) { +Manifold::Impl::Impl(Shape shape, const mat3x4 m) { std::vector vertPos; std::vector triVerts; switch (shape) { @@ -515,9 +515,9 @@ void Manifold::Impl::WarpBatch(std::function)> warpFunc) { Finish(); } -Manifold::Impl Manifold::Impl::Transform(const mat4x3& transform_) const { +Manifold::Impl Manifold::Impl::Transform(const mat3x4& transform_) const { ZoneScoped; - if (transform_ == mat4x3(1.0)) return *this; + if (transform_ == Identity3x4()) return *this; auto policy = autoPolicy(NumVert()); Impl result; if (status_ != Manifold::Error::NoError) { @@ -533,7 +533,7 @@ Manifold::Impl Manifold::Impl::Transform(const mat4x3& transform_) const { result.meshRelation_.originalID = -1; for (auto& m : result.meshRelation_.meshIDtransform) { - m.second.transform = transform_ * mat4(m.second.transform); + m.second.transform = transform_ * Mat4(m.second.transform); } result.vertPos_.resize(NumVert()); @@ -548,7 +548,7 @@ Manifold::Impl Manifold::Impl::Transform(const mat4x3& transform_) const { transform(vertNormal_.begin(), vertNormal_.end(), result.vertNormal_.begin(), TransformNormals({normalTransform})); - const bool invert = glm::determinant(mat3(transform_)) < 0; + const bool invert = la::determinant(mat3(transform_)) < 0; if (halfedgeTangent_.size() > 0) { for_each_n(policy, countAt(0), halfedgeTangent_.size(), @@ -597,7 +597,7 @@ void Manifold::Impl::CalculateNormals() { ZoneScoped; vertNormal_.resize(NumVert()); auto policy = autoPolicy(NumTri(), 1e4); - fill(vertNormal_.begin(), vertNormal_.end(), vec3(0)); + fill(vertNormal_.begin(), vertNormal_.end(), vec3(0.0)); bool calculateTriNormal = false; if (faceNormal_.size() != NumTri()) { faceNormal_.resize(NumTri()); diff --git a/src/impl.h b/src/impl.h index 45e3c95f4..1ba94ad3e 100644 --- a/src/impl.h +++ b/src/impl.h @@ -28,7 +28,7 @@ namespace manifold { struct Manifold::Impl { struct Relation { int originalID = -1; - mat4x3 transform = mat4x3(1); + mat3x4 transform = Identity3x4(); bool backSide = false; }; struct MeshRelationD { @@ -60,7 +60,7 @@ struct Manifold::Impl { Impl() {} enum class Shape { Tetrahedron, Cube, Octahedron }; - Impl(Shape, const mat4x3 = mat4x3(1)); + Impl(Shape, const mat3x4 = Identity3x4()); template Impl(const MeshGLP& meshGL) { @@ -153,10 +153,11 @@ struct Manifold::Impl { meshRelation_.meshIDtransform[meshID] = {originalID}; } else { const Precision* m = meshGL.runTransform.data() + 12 * i; - meshRelation_.meshIDtransform[meshID] = { - originalID, - {m[0], m[1], m[2], m[3], m[4], m[5], m[6], m[7], m[8], m[9], - m[10], m[11]}}; + meshRelation_.meshIDtransform[meshID] = {originalID, + {{m[0], m[1], m[2]}, + {m[3], m[4], m[5]}, + {m[6], m[7], m[8]}, + {m[9], m[10], m[11]}}}; } } } @@ -250,7 +251,7 @@ struct Manifold::Impl { void MarkFailure(Error status); void Warp(std::function warpFunc); void WarpBatch(std::function)> warpFunc); - Impl Transform(const mat4x3& transform) const; + Impl Transform(const mat3x4& transform) const; SparseIndices EdgeCollisions(const Impl& B, bool inverted = false) const; SparseIndices VertexCollisionsZ(VecView vertsIn, bool inverted = false) const; @@ -292,7 +293,7 @@ struct Manifold::Impl { void Face2Tri(const Vec& faceEdge, const Vec& halfedgeRef); PolygonsIdx Face2Polygons(VecView::IterC start, VecView::IterC end, - mat3x2 projection) const; + mat2x3 projection) const; Polygons Slice(double height) const; Polygons Project() const; diff --git a/src/manifold.cpp b/src/manifold.cpp index f29387950..d91ed5b14 100644 --- a/src/manifold.cpp +++ b/src/manifold.cpp @@ -47,13 +47,13 @@ struct UpdateProperties { }; Manifold Halfspace(Box bBox, vec3 normal, double originOffset) { - normal = glm::normalize(normal); + normal = la::normalize(normal); Manifold cutter = Manifold::Cube(vec3(2.0), true).Translate({1.0, 0.0, 0.0}); - double size = glm::length(bBox.Center() - normal * originOffset) + - 0.5 * glm::length(bBox.Size()); + double size = la::length(bBox.Center() - normal * originOffset) + + 0.5 * la::length(bBox.Size()); cutter = cutter.Scale(vec3(size)).Translate({originOffset, 0.0, 0.0}); - double yDeg = glm::degrees(-std::asin(normal.z)); - double zDeg = glm::degrees(std::atan2(normal.y, normal.x)); + double yDeg = degrees(-std::asin(normal.z)); + double zDeg = degrees(std::atan2(normal.y, normal.x)); return cutter.Rotate(0.0, yDeg, zDeg); } @@ -67,7 +67,7 @@ MeshGLP GetMeshGLImpl(const manifold::Manifold::Impl& impl, const bool isOriginal = impl.meshRelation_.originalID >= 0; const bool updateNormals = - !isOriginal && glm::all(glm::greaterThan(normalIdx, ivec3(2))); + !isOriginal && la::all(la::greater(normalIdx, ivec3(2))); MeshGLP out; out.precision = @@ -197,7 +197,7 @@ MeshGLP GetMeshGLImpl(const manifold::Manifold::Impl& impl, for (int i : {0, 1, 2}) { normal[i] = out.vertProperties[start + normalIdx[i]]; } - normal = glm::normalize(runNormalTransform[run] * normal); + normal = la::normalize(runNormalTransform[run] * normal); for (int i : {0, 1, 2}) { out.vertProperties[start + normalIdx[i]] = normal[i]; } @@ -532,7 +532,7 @@ Manifold Manifold::Rotate(double xDegrees, double yDegrees, * * @param m The affine transform matrix to apply to all the vertices. */ -Manifold Manifold::Transform(const mat4x3& m) const { +Manifold Manifold::Transform(const mat3x4& m) const { return Manifold(pNode_->Transform(m)); } @@ -545,11 +545,11 @@ Manifold Manifold::Transform(const mat4x3& m) const { * @param normal The normal vector of the plane to be mirrored over */ Manifold Manifold::Mirror(vec3 normal) const { - if (glm::length(normal) == 0.) { + if (la::length(normal) == 0.) { return Manifold(); } - auto n = glm::normalize(normal); - auto m = mat4x3(mat3(1.0) - 2.0 * glm::outerProduct(n, n)); + auto n = la::normalize(normal); + auto m = mat3x4(mat3(la::identity) - 2.0 * la::outerprod(n, n), vec3()); return Manifold(pNode_->Transform(m)); } @@ -780,7 +780,7 @@ Manifold Manifold::RefineToLength(double length) const { length = std::abs(length); auto pImpl = std::make_shared(*GetCsgLeafNode().GetImpl()); pImpl->Refine([length](vec3 edge, vec4 tangentStart, vec4 tangentEnd) { - return static_cast(glm::length(edge) / length); + return static_cast(la::length(edge) / length); }); return Manifold(std::make_shared(pImpl)); } @@ -803,16 +803,16 @@ Manifold Manifold::RefineToPrecision(double precision) const { if (!pImpl->halfedgeTangent_.empty()) { pImpl->Refine( [precision](vec3 edge, vec4 tangentStart, vec4 tangentEnd) { - const vec3 edgeNorm = glm::normalize(edge); + const vec3 edgeNorm = la::normalize(edge); // Weight heuristic const vec3 tStart = vec3(tangentStart); const vec3 tEnd = vec3(tangentEnd); // Perpendicular to edge - const vec3 start = tStart - edgeNorm * glm::dot(edgeNorm, tStart); - const vec3 end = tEnd - edgeNorm * glm::dot(edgeNorm, tEnd); + const vec3 start = tStart - edgeNorm * la::dot(edgeNorm, tStart); + const vec3 end = tEnd - edgeNorm * la::dot(edgeNorm, tEnd); // Circular arc result plus heuristic term for non-circular curves - const double d = 0.5 * (glm::length(start) + glm::length(end)) + - glm::length(start - end); + const double d = 0.5 * (la::length(start) + la::length(end)) + + la::length(start - end); return static_cast(std::sqrt(3 * d / (4 * precision))); }, true); diff --git a/src/meshIO/meshIO.cpp b/src/meshIO/meshIO.cpp index a985d4b7d..090931416 100644 --- a/src/meshIO/meshIO.cpp +++ b/src/meshIO/meshIO.cpp @@ -45,7 +45,7 @@ aiScene* CreateScene(const ExportOptions& options) { material->AddProperty(&options.mat.roughness, 1, AI_MATKEY_ROUGHNESS_FACTOR); material->AddProperty(&options.mat.metalness, 1, AI_MATKEY_METALLIC_FACTOR); const vec4& color = options.mat.color; - aiColor4D baseColor(color.r, color.g, color.b, color.a); + aiColor4D baseColor(color.x, color.y, color.z, color.w); material->AddProperty(&baseColor, 1, AI_MATKEY_COLOR_DIFFUSE); scene->mNumMeshes = 1; @@ -55,7 +55,7 @@ aiScene* CreateScene(const ExportOptions& options) { scene->mRootNode = new aiNode(); scene->mRootNode->mNumMeshes = 1; - scene->mRootNode->mMeshes = new glm::uint[scene->mRootNode->mNumMeshes]; + scene->mRootNode->mMeshes = new uint[scene->mRootNode->mNumMeshes]; scene->mRootNode->mMeshes[0] = 0; scene->mMeshes[0]->mPrimitiveTypes = aiPrimitiveType_TRIANGLE; @@ -240,8 +240,8 @@ void ExportMesh(const std::string& filename, const MeshGL& mesh, ? 1 : mesh.vertProperties[i * mesh.numProp + options.mat.colorChannels[j]]; - c = glm::saturate(c); - mesh_out->mColors[0][i] = aiColor4D(c.r, c.g, c.b, c.a); + c = la::clamp(c, 0, 1); + mesh_out->mColors[0][i] = aiColor4D(c.x, c.y, c.z, c.w); } } @@ -251,7 +251,7 @@ void ExportMesh(const std::string& filename, const MeshGL& mesh, for (size_t i = 0; i < mesh_out->mNumFaces; ++i) { aiFace& face = mesh_out->mFaces[i]; face.mNumIndices = 3; - face.mIndices = new glm::uint[face.mNumIndices]; + face.mIndices = new uint[face.mNumIndices]; for (int j : {0, 1, 2}) face.mIndices[j] = mesh.triVerts[3 * i + j]; } diff --git a/src/mesh_fixes.h b/src/mesh_fixes.h index 57cb6a650..d4aa800b4 100644 --- a/src/mesh_fixes.h +++ b/src/mesh_fixes.h @@ -27,7 +27,7 @@ struct TransformNormals { mat3 transform; vec3 operator()(vec3 normal) const { - normal = glm::normalize(transform * normal); + normal = la::normalize(transform * normal); if (std::isnan(normal.x)) normal = vec3(0.0); return normal; } diff --git a/src/polygon.cpp b/src/polygon.cpp index 66dc8e8ff..8f8fcdaa9 100644 --- a/src/polygon.cpp +++ b/src/polygon.cpp @@ -30,7 +30,7 @@ static ExecutionParams params; constexpr double kBest = -std::numeric_limits::infinity(); -// it seems that MSVC cannot optimize glm::determinant(mat2(a, b)) +// it seems that MSVC cannot optimize la::determinant(mat2(a, b)) constexpr double determinant2x2(vec2 a, vec2 b) { return a.x * b.y - a.y * b.x; } @@ -176,15 +176,15 @@ bool IsConvex(const PolygonsIdx &polys, double precision) { // Zero-length edges comes out NaN, which won't trip the early return, but // it's okay because that zero-length edge will also get tested // non-normalized and will trip det == 0. - vec2 lastEdge = glm::normalize(firstEdge); + vec2 lastEdge = la::normalize(firstEdge); for (size_t v = 0; v < poly.size(); ++v) { const vec2 edge = v + 1 < poly.size() ? poly[v + 1].pos - poly[v].pos : firstEdge; const double det = determinant2x2(lastEdge, edge); if (det <= 0 || - (std::abs(det) < precision && glm::dot(lastEdge, edge) < 0)) + (std::abs(det) < precision && la::dot(lastEdge, edge) < 0)) return false; - lastEdge = glm::normalize(edge); + lastEdge = la::normalize(edge); } } return true; @@ -326,7 +326,7 @@ class EarClip { // cause CW triangles that exceed precision due to rounding error. bool IsShort(double precision) const { const vec2 edge = right->pos - pos; - return glm::dot(edge, edge) * 4 < precision * precision; + return la::dot(edge, edge) * 4 < precision * precision; } // Like CCW, returns 1 if v is on the inside of the angle formed at this @@ -334,7 +334,7 @@ class EarClip { // Ensure v is more than precision from pos, as case this will not return 0. int Interior(vec2 v, double precision) const { const vec2 diff = v - pos; - if (glm::dot(diff, diff) < precision * precision) { + if (la::dot(diff, diff) < precision * precision) { return 0; } return CCW(pos, left->pos, right->pos, precision) + @@ -357,21 +357,21 @@ class EarClip { while (nextL != nextR && tail != nextR && nextL != (toLeft ? right : left)) { const vec2 edgeL = nextL->pos - center->pos; - const double l2 = glm::dot(edgeL, edgeL); + const double l2 = la::dot(edgeL, edgeL); if (l2 <= p2) { nextL = toLeft ? nextL->left : nextL->right; continue; } const vec2 edgeR = nextR->pos - center->pos; - const double r2 = glm::dot(edgeR, edgeR); + const double r2 = la::dot(edgeR, edgeR); if (r2 <= p2) { nextR = nextR->right; continue; } const vec2 vecLR = nextR->pos - nextL->pos; - const double lr2 = glm::dot(vecLR, vecLR); + const double lr2 = la::dot(vecLR, vecLR); if (lr2 <= p2) { last = center; center = nextL; @@ -411,7 +411,7 @@ class EarClip { if (convexity != 0) { return convexity > 0; } - if (glm::dot(left->pos - pos, right->pos - pos) <= 0) { + if (la::dot(left->pos - pos, right->pos - pos) <= 0) { return true; } return left->InsideEdge(left->right, precision, true); @@ -429,7 +429,7 @@ class EarClip { // right of start - all within a precision tolerance. If onTop != 0, this // restricts which end is allowed to terminate within the precision band. double InterpY2X(vec2 start, int onTop, double precision) const { - if (glm::abs(pos.y - start.y) <= precision) { + if (la::abs(pos.y - start.y) <= precision) { if (right->pos.y <= start.y + precision || onTop == 1) { return NAN; } else { @@ -478,7 +478,7 @@ class EarClip { // to aid in prioritization and produce cleaner triangulations. This doesn't // affect robustness, but may be adjusted to improve output. static double DelaunayCost(vec2 diff, double scale, double precision) { - return -precision - scale * glm::dot(diff, diff); + return -precision - scale * la::dot(diff, diff); } // This is the expensive part of the algorithm, checking this ear against @@ -494,11 +494,11 @@ class EarClip { double EarCost(double precision, IdxCollider &collider) const { vec2 openSide = left->pos - right->pos; const vec2 center = 0.5 * (left->pos + right->pos); - const double scale = 4 / glm::dot(openSide, openSide); - const double radius = glm::length(openSide) / 2; - openSide = glm::normalize(openSide); + const double scale = 4 / la::dot(openSide, openSide); + const double radius = la::length(openSide) / 2; + openSide = la::normalize(openSide); - double totalCost = glm::dot(left->rightDir, rightDir) - 1 - precision; + double totalCost = la::dot(left->rightDir, rightDir) - 1 - precision; if (CCW(pos, left->pos, right->pos, precision) == 0) { // Clip folded ears first return totalCost; @@ -544,7 +544,7 @@ class EarClip { }; static vec2 SafeNormalize(vec2 v) { - vec2 n = glm::normalize(v); + vec2 n = la::normalize(v); return std::isfinite(n.x) ? n : vec2(0, 0); } @@ -616,7 +616,7 @@ class EarClip { } if (ear->IsShort(precision_) || (CCW(ear->left->pos, ear->pos, ear->right->pos, precision_) == 0 && - glm::dot(ear->left->pos - ear->pos, ear->right->pos - ear->pos) > 0 && + la::dot(ear->left->pos - ear->pos, ear->right->pos - ear->pos) > 0 && ear->IsConvex(precision_))) { ClipEar(ear); ClipIfDegenerate(ear->left); @@ -759,7 +759,7 @@ class EarClip { : edge->right->pos.y - start->pos.y > start->pos.y - edge->pos.y ? edge : edge->right; - if (glm::abs(connector->pos.y - start->pos.y) <= precision_) { + if (la::abs(connector->pos.y - start->pos.y) <= precision_) { return connector; } const double above = connector->pos.y > start->pos.y ? 1 : -1; diff --git a/src/properties.cpp b/src/properties.cpp index 36eca8209..ab42c597b 100644 --- a/src/properties.cpp +++ b/src/properties.cpp @@ -33,12 +33,12 @@ struct FaceAreaVolume { const int j = (i + 1) % 3; edge[i] = vertPos[halfedges[3 * face + j].startVert] - vertPos[halfedges[3 * face + i].startVert]; - perimeter += glm::length(edge[i]); + perimeter += la::length(edge[i]); } - vec3 crossP = glm::cross(edge[0], edge[1]); + vec3 crossP = la::cross(edge[0], edge[1]); - double area = glm::length(crossP); - double volume = glm::dot(crossP, vertPos[halfedges[3 * face].startVert]); + double area = la::length(crossP); + double volume = la::dot(crossP, vertPos[halfedges[3 * face].startVert]); return std::make_pair(area / 2.0, volume / 6.0); } @@ -60,24 +60,24 @@ struct CurvatureAngles { const int startVert = halfedge[3 * tri + i].startVert; const int endVert = halfedge[3 * tri + i].endVert; edge[i] = vertPos[endVert] - vertPos[startVert]; - edgeLength[i] = glm::length(edge[i]); + edgeLength[i] = la::length(edge[i]); edge[i] /= edgeLength[i]; const int neighborTri = halfedge[3 * tri + i].pairedHalfedge / 3; const double dihedral = 0.25 * edgeLength[i] * - std::asin(glm::dot(glm::cross(triNormal[tri], triNormal[neighborTri]), - edge[i])); + std::asin(la::dot(la::cross(triNormal[tri], triNormal[neighborTri]), + edge[i])); AtomicAdd(meanCurvature[startVert], dihedral); AtomicAdd(meanCurvature[endVert], dihedral); AtomicAdd(degree[startVert], 1.0); } vec3 phi; - phi[0] = std::acos(-glm::dot(edge[2], edge[0])); - phi[1] = std::acos(-glm::dot(edge[0], edge[1])); - phi[2] = glm::pi() - phi[0] - phi[1]; + phi[0] = std::acos(-la::dot(edge[2], edge[0])); + phi[1] = std::acos(-la::dot(edge[0], edge[1])); + phi[2] = kPi - phi[0] - phi[1]; const double area3 = edgeLength[0] * edgeLength[1] * - glm::length(glm::cross(edge[0], edge[1])) / 6; + la::length(la::cross(edge[0], edge[1])) / 6; for (int i : {0, 1, 2}) { const int vert = halfedge[3 * tri + i].startVert; @@ -155,7 +155,7 @@ struct CheckCCW { bool operator()(size_t face) const { if (halfedges[3 * face].pairedHalfedge < 0) return true; - const mat3x2 projection = GetAxisAlignedProjection(triNormal[face]); + const mat2x3 projection = GetAxisAlignedProjection(triNormal[face]); vec2 v[3]; for (int i : {0, 1, 2}) v[i] = projection * vertPos[halfedges[3 * face + i].startVert]; @@ -168,12 +168,12 @@ struct CheckCCW { vec2 v1 = v[1] - v[0]; vec2 v2 = v[2] - v[0]; double area = v1.x * v2.y - v1.y * v2.x; - double base2 = std::max(glm::dot(v1, v1), glm::dot(v2, v2)); + double base2 = std::max(la::dot(v1, v1), la::dot(v2, v2)); double base = std::sqrt(base2); vec3 V0 = vertPos[halfedges[3 * face].startVert]; vec3 V1 = vertPos[halfedges[3 * face + 1].startVert]; vec3 V2 = vertPos[halfedges[3 * face + 2].startVert]; - vec3 norm = glm::cross(V1 - V0, V2 - V0); + vec3 norm = la::cross(V1 - V0, V2 - V0); printf( "Tri %ld does not match normal, approx height = %g, base = %g\n" "tol = %g, area2 = %g, base2*tol2 = %g\n" @@ -271,7 +271,7 @@ void Manifold::Impl::CalculateCurvature(int gaussianIdx, int meanIdx) { if (IsEmpty()) return; if (gaussianIdx < 0 && meanIdx < 0) return; Vec vertMeanCurvature(NumVert(), 0); - Vec vertGaussianCurvature(NumVert(), glm::two_pi()); + Vec vertGaussianCurvature(NumVert(), kTwoPi); Vec vertArea(NumVert(), 0); Vec degree(NumVert(), 0); auto policy = autoPolicy(NumTri(), 1e4); @@ -314,14 +314,14 @@ void Manifold::Impl::CalculateBBox() { vec3(std::numeric_limits::infinity()), [](auto a, auto b) { if (std::isnan(a.x)) return b; if (std::isnan(b.x)) return a; - return glm::min(a, b); + return la::min(a, b); }); bBox_.max = reduce(vertPos_.begin(), vertPos_.end(), vec3(-std::numeric_limits::infinity()), [](auto a, auto b) { if (std::isnan(a.x)) return b; if (std::isnan(b.x)) return a; - return glm::max(a, b); + return la::max(a, b); }); } @@ -333,7 +333,7 @@ bool Manifold::Impl::IsFinite() const { return transform_reduce( vertPos_.begin(), vertPos_.end(), true, [](bool a, bool b) { return a && b; }, - [](auto v) { return glm::all(glm::isfinite(v)); }); + [](auto v) { return la::all(la::isfinite(v)); }); } /** diff --git a/src/quickhull.cpp b/src/quickhull.cpp index 8cb2c577d..14609ce78 100644 --- a/src/quickhull.cpp +++ b/src/quickhull.cpp @@ -29,17 +29,17 @@ double defaultEps() { return 0.0000001; } inline double getSquaredDistanceBetweenPointAndRay(const vec3& p, const Ray& r) { const vec3 s = p - r.S; - double t = glm::dot(s, r.V); - return glm::dot(s, s) - t * t * r.VInvLengthSquared; + double t = la::dot(s, r.V); + return la::dot(s, s) - t * t * r.VInvLengthSquared; } inline double getSquaredDistance(const vec3& p1, const vec3& p2) { - return glm::dot(p1 - p2, p1 - p2); + return la::dot(p1 - p2, p1 - p2); } // Note that the unit of distance returned is relative to plane's normal's // length (divide by N.getNormalized() if needed to get the "real" distance). inline double getSignedDistanceToPlane(const vec3& v, const Plane& p) { - return glm::dot(p.N, v) + p.D; + return la::dot(p.N, v) + p.D; } inline vec3 getTriangleNormal(const vec3& a, const vec3& b, const vec3& c) { @@ -53,7 +53,7 @@ inline vec3 getTriangleNormal(const vec3& a, const vec3& b, const vec3& c) { double px = y * rhsz - z * rhsy; double py = z * rhsx - x * rhsz; double pz = x * rhsy - y * rhsx; - return glm::normalize(vec3(px, py, pz)); + return la::normalize(vec3(px, py, pz)); } size_t MeshBuilder::addFace() { @@ -385,7 +385,7 @@ void QuickHull::createConvexHalfedgeMesh() { } else { const Plane& P = pvf.P; pvf.visibilityCheckedOnIteration = iter; - const double d = glm::dot(P.N, activePoint) + P.D; + const double d = la::dot(P.N, activePoint) + P.D; if (d > 0) { pvf.isVisibleFaceOnCurrentIteration = 1; pvf.horizonEdgesOnCurrentIteration = 0; diff --git a/src/quickhull.h b/src/quickhull.h index 729de29b2..3f9057509 100644 --- a/src/quickhull.h +++ b/src/quickhull.h @@ -95,7 +95,7 @@ class Plane { double sqrNLength; bool isPointOnPositiveSide(const vec3& Q) const { - double d = glm::dot(N, Q) + D; + double d = la::dot(N, Q) + D; if (d >= 0) return true; return false; } @@ -104,7 +104,7 @@ class Plane { // Construct a plane using normal N and any point P on the plane Plane(const vec3& N, const vec3& P) - : N(N), D(glm::dot(-N, P)), sqrNLength(glm::dot(N, N)) {} + : N(N), D(la::dot(-N, P)), sqrNLength(la::dot(N, N)) {} }; struct Ray { @@ -113,7 +113,7 @@ struct Ray { const double VInvLengthSquared; Ray(const vec3& S, const vec3& V) - : S(S), V(V), VInvLengthSquared(1 / (glm::dot(V, V))) {} + : S(S), V(V), VInvLengthSquared(1 / (la::dot(V, V))) {} }; class MeshBuilder { diff --git a/src/sdf.cpp b/src/sdf.cpp index fb3ef1ce5..ab806025a 100644 --- a/src/sdf.cpp +++ b/src/sdf.cpp @@ -12,8 +12,6 @@ // See the License for the specific language governing permissions and // limitations under the License. -#include - #include "./hashtable.h" #include "./impl.h" #include "./utils.h" @@ -148,15 +146,15 @@ inline vec3 FindSurface(vec3 pos0, double d0, vec3 pos1, double d1, double tol, // Sole tuning parameter, k: (0, 1) - smaller value gets better median // performance, but also hits the worst case more often. const double k = 0.1; - const double check = 2 * tol / glm::length(pos0 - pos1); + const double check = 2 * tol / la::length(pos0 - pos1); double frac = 1; double biFrac = 1; while (frac > check) { - const double t = glm::mix(d0 / (d0 - d1), 0.5, k); + const double t = la::lerp(d0 / (d0 - d1), 0.5, k); const double r = biFrac / frac - 0.5; - const double x = glm::abs(t - 0.5) < r ? t : 0.5 - r * (t < 0.5 ? 1 : -1); + const double x = la::abs(t - 0.5) < r ? t : 0.5 - r * (t < 0.5 ? 1 : -1); - const vec3 mid = glm::mix(pos0, pos1, x); + const vec3 mid = la::lerp(pos0, pos1, x); const double d = sdf(mid) - level; if ((d > 0) == (d0 > 0)) { @@ -171,7 +169,7 @@ inline vec3 FindSurface(vec3 pos0, double d0, vec3 pos1, double d1, double tol, biFrac /= 2; } - return glm::mix(pos0, pos1, d0 / (d0 - d1)); + return la::lerp(pos0, pos1, d0 / (d0 - d1)); } /** @@ -217,7 +215,7 @@ struct NearSurface { const ivec4 gridIndex = DecodeIndex(index, gridPow); - if (glm::any(glm::greaterThan(ivec3(gridIndex), gridSize))) return; + if (la::any(la::greater(ivec3(gridIndex), gridSize))) return; GridVert gridVert; gridVert.distance = voxels[EncodeIndex(gridIndex + kVoxelOffset, gridPow)]; @@ -239,14 +237,14 @@ struct NearSurface { ++opposedVerts; } // Approximate bound on vert movement. - if (glm::abs(val) > kD * glm::abs(gridVert.distance) && - glm::abs(val) > glm::abs(vMax)) { + if (la::abs(val) > kD * la::abs(gridVert.distance) && + la::abs(val) > la::abs(vMax)) { vMax = val; closestNeighbor = i; } } else if (!gridVert.SameSide(valOp) && - glm::abs(valOp) > kD * glm::abs(gridVert.distance) && - glm::abs(valOp) > glm::abs(vMax)) { + la::abs(valOp) > kD * la::abs(gridVert.distance) && + la::abs(valOp) > la::abs(vMax)) { vMax = valOp; closestNeighbor = i + 7; } @@ -264,7 +262,7 @@ struct NearSurface { Position(neighborIndex, origin, spacing), vMax, tol, level, sdf); // Bound the delta of each vert to ensure the tetrahedron cannot invert. - if (glm::all(glm::lessThan(glm::abs(pos - gridPos), kS * spacing))) { + if (la::all(la::less(la::abs(pos - gridPos), kS * spacing))) { const int idx = AtomicAdd(vertIndex[0], 1); vertPos[idx] = pos; gridVert.movedVert = idx; @@ -469,7 +467,7 @@ Manifold Manifold::LevelSet(std::function sdf, Box bounds, const ivec3 gridSize(dim / edgeLength + 1.0); const vec3 spacing = dim / (vec3(gridSize - 1)); - const ivec3 gridPow(glm::log2(gridSize + 2) + 1); + const ivec3 gridPow(la::log2(gridSize + 2) + 1); const Uint64 maxIndex = EncodeIndex(ivec4(gridSize + 2, 1), gridPow); // Parallel policies violate will crash language runtimes with runtime locks @@ -488,7 +486,7 @@ Manifold Manifold::LevelSet(std::function sdf, Box bounds, }); size_t tableSize = std::min( - 2 * maxIndex, static_cast(10 * glm::pow(maxIndex, 0.667))); + 2 * maxIndex, static_cast(10 * la::pow(maxIndex, 0.667))); HashTable gridVerts(tableSize); vertPos.resize(gridVerts.Size() * 7); @@ -501,7 +499,7 @@ Manifold Manifold::LevelSet(std::function sdf, Box bounds, if (gridVerts.Full()) { // Resize HashTable const vec3 lastVert = vertPos[index[0] - 1]; const Uint64 lastIndex = - EncodeIndex(ivec4((lastVert - origin) / spacing, 1), gridPow); + EncodeIndex(ivec4(ivec3((lastVert - origin) / spacing), 1), gridPow); const double ratio = static_cast(maxIndex) / lastIndex; if (ratio > 1000) // do not trust the ratio if it is too large diff --git a/src/shared.h b/src/shared.h index 51f17b93a..83759e0b7 100644 --- a/src/shared.h +++ b/src/shared.h @@ -25,8 +25,8 @@ namespace manifold { * @{ */ inline vec3 SafeNormalize(vec3 v) { - v = glm::normalize(v); - return std::isfinite(v.x) ? v : vec3(0); + v = la::normalize(v); + return std::isfinite(v.x) ? v : vec3(0.0); } inline double MaxPrecision(double minPrecision, const Box& bBox) { @@ -40,52 +40,52 @@ inline int NextHalfedge(int current) { return current; } -inline mat3 NormalTransform(const mat4x3& transform) { - return glm::inverse(glm::transpose(mat3(transform))); +inline mat3 NormalTransform(const mat3x4& transform) { + return la::inverse(la::transpose(mat3(transform))); } /** * By using the closest axis-aligned projection to the normal instead of a * projection along the normal, we avoid introducing any rounding error. */ -inline mat3x2 GetAxisAlignedProjection(vec3 normal) { - vec3 absNormal = glm::abs(normal); +inline mat2x3 GetAxisAlignedProjection(vec3 normal) { + vec3 absNormal = la::abs(normal); double xyzMax; - mat2x3 projection; + mat3x2 projection; if (absNormal.z > absNormal.x && absNormal.z > absNormal.y) { - projection = mat2x3(1.0, 0.0, 0.0, // - 0.0, 1.0, 0.0); + projection = mat3x2({1.0, 0.0, 0.0}, // + {0.0, 1.0, 0.0}); xyzMax = normal.z; } else if (absNormal.y > absNormal.x) { - projection = mat2x3(0.0, 0.0, 1.0, // - 1.0, 0.0, 0.0); + projection = mat3x2({0.0, 0.0, 1.0}, // + {1.0, 0.0, 0.0}); xyzMax = normal.y; } else { - projection = mat2x3(0.0, 1.0, 0.0, // - 0.0, 0.0, 1.0); + projection = mat3x2({0.0, 1.0, 0.0}, // + {0.0, 0.0, 1.0}); xyzMax = normal.x; } if (xyzMax < 0) projection[0] *= -1.0; - return glm::transpose(projection); + return la::transpose(projection); } inline vec3 GetBarycentric(const vec3& v, const mat3& triPos, double precision) { const mat3 edges(triPos[2] - triPos[1], triPos[0] - triPos[2], triPos[1] - triPos[0]); - const vec3 d2(glm::dot(edges[0], edges[0]), glm::dot(edges[1], edges[1]), - glm::dot(edges[2], edges[2])); + const vec3 d2(la::dot(edges[0], edges[0]), la::dot(edges[1], edges[1]), + la::dot(edges[2], edges[2])); const int longSide = d2[0] > d2[1] && d2[0] > d2[2] ? 0 : d2[1] > d2[2] ? 1 : 2; - const vec3 crossP = glm::cross(edges[0], edges[1]); - const double area2 = glm::dot(crossP, crossP); + const vec3 crossP = la::cross(edges[0], edges[1]); + const double area2 = la::dot(crossP, crossP); const double tol2 = precision * precision; - vec3 uvw(0); + vec3 uvw(0.0); for (const int i : {0, 1, 2}) { const vec3 dv = v - triPos[i]; - if (glm::dot(dv, dv) < tol2) { + if (la::dot(dv, dv) < tol2) { // Return exactly equal if within tolerance of vert. uvw[i] = 1; return uvw; @@ -97,17 +97,17 @@ inline vec3 GetBarycentric(const vec3& v, const mat3& triPos, } else if (area2 > d2[longSide] * tol2) { // triangle for (const int i : {0, 1, 2}) { const int j = Next3(i); - const vec3 crossPv = glm::cross(edges[i], v - triPos[j]); - const double area2v = glm::dot(crossPv, crossPv); + const vec3 crossPv = la::cross(edges[i], v - triPos[j]); + const double area2v = la::dot(crossPv, crossPv); // Return exactly equal if within tolerance of edge. - uvw[i] = area2v < d2[i] * tol2 ? 0 : glm::dot(crossPv, crossP); + uvw[i] = area2v < d2[i] * tol2 ? 0 : la::dot(crossPv, crossP); } uvw /= (uvw[0] + uvw[1] + uvw[2]); return uvw; } else { // line const int nextV = Next3(longSide); const double alpha = - glm::dot(v - triPos[nextV], edges[longSide]) / d2[longSide]; + la::dot(v - triPos[nextV], edges[longSide]) / d2[longSide]; uvw[longSide] = 0; uvw[nextV] = 1 - alpha; const int lastV = Next3(nextV); diff --git a/src/smoothing.cpp b/src/smoothing.cpp index b8d28fde1..c0e4aaeab 100644 --- a/src/smoothing.cpp +++ b/src/smoothing.cpp @@ -12,9 +12,6 @@ // See the License for the specific language governing permissions and // limitations under the License. -#define GLM_ENABLE_EXPERIMENTAL -#include - #include "./impl.h" #include "manifold/parallel.h" @@ -25,23 +22,23 @@ using namespace manifold; // unless in and ref are colinear, in which case it falls back to the plane of // ref and altIn. vec3 OrthogonalTo(vec3 in, vec3 altIn, vec3 ref) { - vec3 out = in - glm::dot(in, ref) * ref; - if (glm::dot(out, out) < kTolerance * glm::dot(in, in)) { - out = altIn - glm::dot(altIn, ref) * ref; + vec3 out = in - la::dot(in, ref) * ref; + if (la::dot(out, out) < kTolerance * la::dot(in, in)) { + out = altIn - la::dot(altIn, ref) * ref; } return SafeNormalize(out); } double Wrap(double radians) { - return radians < -glm::pi() ? radians + glm::two_pi() - : radians > glm::pi() ? radians - glm::two_pi() - : radians; + return radians < -kPi ? radians + kTwoPi + : radians > kPi ? radians - kTwoPi + : radians; } // Get the angle between two unit-vectors. double AngleBetween(vec3 a, vec3 b) { - const double dot = glm::dot(a, b); - return dot >= 1 ? 0 : (dot <= -1 ? glm::pi() : glm::acos(dot)); + const double dot = la::dot(a, b); + return dot >= 1 ? 0 : (dot <= -1 ? kPi : la::acos(dot)); } // Calculate a tangent vector in the form of a weighted cubic Bezier taking as @@ -52,11 +49,11 @@ double AngleBetween(vec3 a, vec3 b) { vec4 CircularTangent(const vec3& tangent, const vec3& edgeVec) { const vec3 dir = SafeNormalize(tangent); - double weight = std::max(0.5, glm::dot(dir, SafeNormalize(edgeVec))); + double weight = std::max(0.5, la::dot(dir, SafeNormalize(edgeVec))); // Quadratic weighted bezier for circular interpolation - const vec4 bz2 = vec4(dir * 0.5 * glm::length(edgeVec), weight); + const vec4 bz2 = vec4(dir * 0.5 * la::length(edgeVec), weight); // Equivalent cubic weighted bezier - const vec4 bz3 = glm::mix(vec4(0, 0, 0, 1), bz2, 2 / 3.0); + const vec4 bz3 = la::lerp(vec4(0, 0, 0, 1), bz2, 2 / 3.0); // Convert from homogeneous form to geometric form return vec4(vec3(bz3) / bz3.w, bz3.w); } @@ -85,30 +82,30 @@ struct InterpTri { return Homogeneous(vec4(point, 0) + tangent); } - static mat2x4 CubicBezier2Linear(vec4 p0, vec4 p1, vec4 p2, vec4 p3, + static mat4x2 CubicBezier2Linear(vec4 p0, vec4 p1, vec4 p2, vec4 p3, double x) { - mat2x4 out; - vec4 p12 = glm::mix(p1, p2, x); - out[0] = glm::mix(glm::mix(p0, p1, x), p12, x); - out[1] = glm::mix(p12, glm::mix(p2, p3, x), x); + mat4x2 out; + vec4 p12 = la::lerp(p1, p2, x); + out[0] = la::lerp(la::lerp(p0, p1, x), p12, x); + out[1] = la::lerp(p12, la::lerp(p2, p3, x), x); return out; } - static vec3 BezierPoint(mat2x4 points, double x) { - return HNormalize(glm::mix(points[0], points[1], x)); + static vec3 BezierPoint(mat4x2 points, double x) { + return HNormalize(la::lerp(points[0], points[1], x)); } - static vec3 BezierTangent(mat2x4 points) { + static vec3 BezierTangent(mat4x2 points) { return SafeNormalize(HNormalize(points[1]) - HNormalize(points[0])); } static vec3 RotateFromTo(vec3 v, quat start, quat end) { - return end * glm::conjugate(start) * v; + return la::qrot(end, la::qrot(la::qconj(start), v)); } static quat Slerp(const quat& x, const quat& y, double a, bool longWay) { quat z = y; - double cosTheta = glm::dot(x, y); + double cosTheta = la::dot(x, y); // Take the long way around the sphere only when requested if ((cosTheta < 0) != longWay) { @@ -116,8 +113,8 @@ struct InterpTri { cosTheta = -cosTheta; } - if (cosTheta > 1.0 - glm::epsilon()) { - return glm::lerp(x, z, a); // for numerical stability + if (cosTheta > 1.0 - std::numeric_limits::epsilon()) { + return la::lerp(x, z, a); // for numerical stability } else { double angle = std::acos(cosTheta); return (std::sin((1.0 - a) * angle) * x + std::sin(a * angle) * z) / @@ -125,51 +122,49 @@ struct InterpTri { } } - static mat2x4 Bezier2Bezier(const mat2x3& corners, const mat2x4& tangentsX, - const mat2x4& tangentsY, double x, + static mat4x2 Bezier2Bezier(const mat3x2& corners, const mat4x2& tangentsX, + const mat4x2& tangentsY, double x, const vec3& anchor) { - const mat2x4 bez = CubicBezier2Linear( + const mat4x2 bez = CubicBezier2Linear( Homogeneous(corners[0]), Bezier(corners[0], tangentsX[0]), Bezier(corners[1], tangentsX[1]), Homogeneous(corners[1]), x); const vec3 end = BezierPoint(bez, x); const vec3 tangent = BezierTangent(bez); - const mat2x3 nTangentsX(SafeNormalize(vec3(tangentsX[0])), + const mat3x2 nTangentsX(SafeNormalize(vec3(tangentsX[0])), -SafeNormalize(vec3(tangentsX[1]))); - const mat2x3 biTangents = { + const mat3x2 biTangents = { OrthogonalTo(vec3(tangentsY[0]), (anchor - corners[0]), nTangentsX[0]), OrthogonalTo(vec3(tangentsY[1]), (anchor - corners[1]), nTangentsX[1])}; - const quat q0 = - glm::quat_cast(mat3(nTangentsX[0], biTangents[0], - glm::cross(nTangentsX[0], biTangents[0]))); - const quat q1 = - glm::quat_cast(mat3(nTangentsX[1], biTangents[1], - glm::cross(nTangentsX[1], biTangents[1]))); + const quat q0 = la::rotation_quat(mat3( + nTangentsX[0], biTangents[0], la::cross(nTangentsX[0], biTangents[0]))); + const quat q1 = la::rotation_quat(mat3( + nTangentsX[1], biTangents[1], la::cross(nTangentsX[1], biTangents[1]))); const vec3 edge = corners[1] - corners[0]; const bool longWay = - glm::dot(nTangentsX[0], edge) + glm::dot(nTangentsX[1], edge) < 0; + la::dot(nTangentsX[0], edge) + la::dot(nTangentsX[1], edge) < 0; const quat qTmp = Slerp(q0, q1, x, longWay); - const quat q = glm::rotation(qTmp * vec3(1, 0, 0), tangent) * qTmp; + const quat q = la::qmul(la::rotation_quat(la::qxdir(qTmp), tangent), qTmp); - const vec3 delta = glm::mix(RotateFromTo(vec3(tangentsY[0]), q0, q), + const vec3 delta = la::lerp(RotateFromTo(vec3(tangentsY[0]), q0, q), RotateFromTo(vec3(tangentsY[1]), q1, q), x); - const double deltaW = glm::mix(tangentsY[0].w, tangentsY[1].w, x); + const double deltaW = la::lerp(tangentsY[0].w, tangentsY[1].w, x); return {Homogeneous(end), vec4(delta, deltaW)}; } - static vec3 Bezier2D(const mat4x3& corners, const mat4& tangentsX, + static vec3 Bezier2D(const mat3x4& corners, const mat4& tangentsX, const mat4& tangentsY, double x, double y, const vec3& centroid) { - mat2x4 bez0 = + mat4x2 bez0 = Bezier2Bezier({corners[0], corners[1]}, {tangentsX[0], tangentsX[1]}, {tangentsY[0], tangentsY[1]}, x, centroid); - mat2x4 bez1 = + mat4x2 bez1 = Bezier2Bezier({corners[2], corners[3]}, {tangentsX[2], tangentsX[3]}, {tangentsY[2], tangentsY[3]}, 1 - x, centroid); - const mat2x4 bez = + const mat4x2 bez = CubicBezier2Linear(bez0[0], Bezier(vec3(bez0[0]), bez0[1]), Bezier(vec3(bez1[0]), bez1[1]), bez1[0], y); return BezierPoint(bez, y); @@ -181,12 +176,12 @@ struct InterpTri { const vec4 uvw = vertBary[vert].uvw; const ivec4 halfedges = impl->GetHalfedges(tri); - const mat4x3 corners = { + const mat3x4 corners = { impl->vertPos_[impl->halfedge_[halfedges[0]].startVert], impl->vertPos_[impl->halfedge_[halfedges[1]].startVert], impl->vertPos_[impl->halfedge_[halfedges[2]].startVert], halfedges[3] < 0 - ? vec3(0) + ? vec3(0.0) : impl->vertPos_[impl->halfedge_[halfedges[3]].startVert]}; for (const int i : {0, 1, 2, 3}) { @@ -196,13 +191,13 @@ struct InterpTri { } } - vec4 posH(0); + vec4 posH(0.0); if (halfedges[3] < 0) { // tri - const mat3x4 tangentR = {impl->halfedgeTangent_[halfedges[0]], + const mat4x3 tangentR = {impl->halfedgeTangent_[halfedges[0]], impl->halfedgeTangent_[halfedges[1]], impl->halfedgeTangent_[halfedges[2]]}; - const mat3x4 tangentL = { + const mat4x3 tangentL = { impl->halfedgeTangent_[impl->halfedge_[halfedges[2]].pairedHalfedge], impl->halfedgeTangent_[impl->halfedge_[halfedges[0]].pairedHalfedge], impl->halfedgeTangent_[impl->halfedge_[halfedges[1]].pairedHalfedge]}; @@ -213,13 +208,13 @@ struct InterpTri { const int k = Prev3(i); const double x = uvw[k] / (1 - uvw[i]); - const mat2x4 bez = + const mat4x2 bez = Bezier2Bezier({corners[j], corners[k]}, {tangentR[j], tangentL[k]}, {tangentL[j], tangentR[k]}, x, centroid); - const mat2x4 bez1 = CubicBezier2Linear( + const mat4x2 bez1 = CubicBezier2Linear( bez[0], Bezier(vec3(bez[0]), bez[1]), - Bezier(corners[i], glm::mix(tangentR[i], tangentL[i], x)), + Bezier(corners[i], la::lerp(tangentR[i], tangentL[i], x)), Homogeneous(corners[i]), uvw[i]); const vec3 p = BezierPoint(bez1, uvw[i]); posH += Homogeneous(vec4(p, uvw[j] * uvw[k])); @@ -279,7 +274,7 @@ vec4 Manifold::Impl::TangentFromNormal(const vec3& normal, int halfedge) const { const vec3 edgeVec = vertPos_[edge.endVert] - vertPos_[edge.startVert]; const vec3 edgeNormal = faceNormal_[halfedge / 3] + faceNormal_[edge.pairedHalfedge / 3]; - vec3 dir = glm::cross(glm::cross(edgeNormal, edgeVec), normal); + vec3 dir = la::cross(la::cross(edgeNormal, edgeVec), normal); return CircularTangent(dir, edgeVec); } @@ -403,12 +398,12 @@ Vec Manifold::Impl::VertHalfedge() const { std::vector Manifold::Impl::SharpenEdges( double minSharpAngle, double minSmoothness) const { std::vector sharpenedEdges; - const double minRadians = glm::radians(minSharpAngle); + const double minRadians = radians(minSharpAngle); for (size_t e = 0; e < halfedge_.size(); ++e) { if (!halfedge_[e].IsForward()) continue; const size_t pair = halfedge_[e].pairedHalfedge; const double dihedral = - std::acos(glm::dot(faceNormal_[e / 3], faceNormal_[pair / 3])); + std::acos(la::dot(faceNormal_[e / 3], faceNormal_[pair / 3])); if (dihedral > minRadians) { sharpenedEdges.push_back({e, minSmoothness}); sharpenedEdges.push_back({pair, minSmoothness}); @@ -449,7 +444,7 @@ void Manifold::Impl::SetNormals(int normalIdx, double minSharpAngle) { const int tri1 = e / 3; const int tri2 = pair / 3; const double dihedral = - glm::degrees(std::acos(glm::dot(faceNormal_[tri1], faceNormal_[tri2]))); + degrees(std::acos(la::dot(faceNormal_[tri1], faceNormal_[tri2]))); if (dihedral > minSharpAngle) { ++vertNumSharp[halfedge_[e].startVert]; ++vertNumSharp[halfedge_[e].endVert]; @@ -520,8 +515,8 @@ void Manifold::Impl::SetNormals(int normalIdx, double minSharpAngle) { int next = NextHalfedge(halfedge_[current].pairedHalfedge); const int face = next / 3; - const double dihedral = glm::degrees( - std::acos(glm::dot(faceNormal_[face], faceNormal_[prevFace]))); + const double dihedral = degrees( + std::acos(la::dot(faceNormal_[face], faceNormal_[prevFace]))); if (dihedral > minSharpAngle || triIsFlatFace[face] != triIsFlatFace[prevFace] || (triIsFlatFace[face] && triIsFlatFace[prevFace] && @@ -566,19 +561,19 @@ void Manifold::Impl::SetNormals(int normalIdx, double minSharpAngle) { }, [this, &triIsFlatFace, &normals, &group, minSharpAngle]( int current, const FaceEdge& here, FaceEdge& next) { - const double dihedral = glm::degrees(std::acos( - glm::dot(faceNormal_[here.face], faceNormal_[next.face]))); + const double dihedral = degrees(std::acos( + la::dot(faceNormal_[here.face], faceNormal_[next.face]))); if (dihedral > minSharpAngle || triIsFlatFace[here.face] != triIsFlatFace[next.face] || (triIsFlatFace[here.face] && triIsFlatFace[next.face] && !meshRelation_.triRef[here.face].SameFace( meshRelation_.triRef[next.face]))) { - normals.push_back(vec3(0)); + normals.push_back(vec3(0.0)); } group.push_back(normals.size() - 1); if (std::isfinite(next.edgeVec.x)) { normals.back() += - SafeNormalize(glm::cross(next.edgeVec, here.edgeVec)) * + SafeNormalize(la::cross(next.edgeVec, here.edgeVec)) * AngleBetween(here.edgeVec, next.edgeVec); } else { next.edgeVec = here.edgeVec; @@ -685,7 +680,7 @@ void Manifold::Impl::DistributeTangents(const Vec& fixedHalfedges) { halfedge = NextHalfedge(halfedge_[halfedge].pairedHalfedge); } - vec3 normal(0); + vec3 normal(0.0); Vec currentAngle; Vec desiredAngle; @@ -704,25 +699,25 @@ void Manifold::Impl::DistributeTangents(const Vec& fixedHalfedges) { SafeNormalize(vertPos_[halfedge_[current].endVert] - center); const vec3 thisTangent = SafeNormalize(vec3(halfedgeTangent_[current])); - normal += glm::cross(thisTangent, lastTangent); + normal += la::cross(thisTangent, lastTangent); // cumulative sum desiredAngle.push_back( AngleBetween(thisEdgeVec, lastEdgeVec) + (desiredAngle.size() > 0 ? desiredAngle.back() : 0)); if (current == halfedge) { - currentAngle.push_back(glm::two_pi()); + currentAngle.push_back(kTwoPi); } else { currentAngle.push_back(AngleBetween(thisTangent, firstTangent)); - if (glm::dot(approxNormal, glm::cross(thisTangent, firstTangent)) < + if (la::dot(approxNormal, la::cross(thisTangent, firstTangent)) < 0) { - currentAngle.back() = glm::two_pi() - currentAngle.back(); + currentAngle.back() = kTwoPi - currentAngle.back(); } } lastEdgeVec = thisEdgeVec; lastTangent = thisTangent; } while (!fixedHalfedges[current]); - if (currentAngle.size() == 1 || glm::dot(normal, normal) == 0) return; + if (currentAngle.size() == 1 || la::dot(normal, normal) == 0) return; const double scale = currentAngle.back() / desiredAngle.back(); double offset = 0; @@ -741,17 +736,17 @@ void Manifold::Impl::DistributeTangents(const Vec& fixedHalfedges) { desiredAngle[i] *= scale; const double lastAngle = i > 0 ? desiredAngle[i - 1] : 0; // shrink obtuse angles - if (desiredAngle[i] - lastAngle > glm::pi()) { - desiredAngle[i] = lastAngle + glm::pi(); + if (desiredAngle[i] - lastAngle > kPi) { + desiredAngle[i] = lastAngle + kPi; } else if (i + 1 < desiredAngle.size() && - scale * desiredAngle[i + 1] - desiredAngle[i] > - glm::pi()) { - desiredAngle[i] = scale * desiredAngle[i + 1] - glm::pi(); + scale * desiredAngle[i + 1] - desiredAngle[i] > kPi) { + desiredAngle[i] = scale * desiredAngle[i + 1] - kPi; } const double angle = currentAngle[i] - desiredAngle[i] - offset; vec3 tangent(halfedgeTangent_[current]); - halfedgeTangent_[current] = vec4(glm::rotate(tangent, angle, normal), - halfedgeTangent_[current].w); + const quat q = la::rotation_quat(la::normalize(normal), angle); + halfedgeTangent_[current] = + vec4(la::qrot(q, tangent), halfedgeTangent_[current].w); ++i; } while (!fixedHalfedges[current]); }); @@ -790,7 +785,7 @@ void Manifold::Impl::CreateTangents(int normalIdx) { const vec3 normal = GetNormal(halfedge, normalIdx); const vec3 diff = faceNormal_[halfedge / 3] - normal; return FlatNormal( - {glm::dot(diff, diff) < kTolerance * kTolerance, normal}); + {la::dot(diff, diff) < kTolerance * kTolerance, normal}); }, [&faceEdges, &tangent, &fixedHalfedge, this]( int halfedge, const FlatNormal& here, const FlatNormal& next) { @@ -801,7 +796,7 @@ void Manifold::Impl::CreateTangents(int normalIdx) { // mark special edges const vec3 diff = next.normal - here.normal; const bool differentNormals = - glm::dot(diff, diff) > kTolerance * kTolerance; + la::dot(diff, diff) > kTolerance * kTolerance; if (differentNormals || here.isFlatFace != next.isFlatFace) { fixedHalfedge[halfedge] = true; if (faceEdges[0] == -1) { @@ -816,9 +811,9 @@ void Manifold::Impl::CreateTangents(int normalIdx) { if (differentNormals) { const vec3 edgeVec = vertPos_[halfedge_[halfedge].endVert] - vertPos_[halfedge_[halfedge].startVert]; - const vec3 dir = glm::cross(here.normal, next.normal); + const vec3 dir = la::cross(here.normal, next.normal); tangent[halfedge] = CircularTangent( - glm::sign(glm::dot(dir, edgeVec)) * dir, edgeVec); + (la::dot(dir, edgeVec) < 0 ? -1.0 : 1.0) * dir, edgeVec); } else { tangent[halfedge] = TangentFromNormal(here.normal, halfedge); } @@ -829,7 +824,7 @@ void Manifold::Impl::CreateTangents(int normalIdx) { vertPos_[halfedge_[faceEdges[0]].startVert]; const vec3 edge1 = vertPos_[halfedge_[faceEdges[1]].endVert] - vertPos_[halfedge_[faceEdges[1]].startVert]; - const vec3 newTangent = glm::normalize(edge0) - glm::normalize(edge1); + const vec3 newTangent = la::normalize(edge0) - la::normalize(edge1); tangent[faceEdges[0]] = CircularTangent(newTangent, edge0); tangent[faceEdges[1]] = CircularTangent(-newTangent, edge1); } else if (faceEdges[0] == -1 && faceEdges[0] == -1) { @@ -933,8 +928,8 @@ void Manifold::Impl::CreateTangents(std::vector sharpenedEdges) { const int second = vert[1].first.halfedge; fixedHalfedge[first] = true; fixedHalfedge[second] = true; - const vec3 newTangent = glm::normalize( - vec3(halfedgeTangent_[first]) - vec3(halfedgeTangent_[second])); + const vec3 newTangent = la::normalize(vec3(halfedgeTangent_[first]) - + vec3(halfedgeTangent_[second])); const vec3 pos = vertPos_[halfedge_[first].startVert]; halfedgeTangent_[first] = CircularTangent( diff --git a/src/subdivision.cpp b/src/subdivision.cpp index 8b1996efd..421e11ada 100644 --- a/src/subdivision.cpp +++ b/src/subdivision.cpp @@ -89,7 +89,7 @@ class Partition { return partition; } - Vec Reindex(ivec4 triVerts, ivec4 edgeOffsets, glm::bvec4 edgeFwd, + Vec Reindex(ivec4 triVerts, ivec4 edgeOffsets, bvec4 edgeFwd, int interiorOffset) const { Vec newVerts; newVerts.reserve(vertBary.size()); @@ -97,7 +97,7 @@ class Partition { ivec4 outTri = {0, 1, 2, 3}; if (triVerts[3] < 0 && idx[1] != Next3(idx[0])) { triIdx = {idx[2], idx[0], idx[1], idx[3]}; - edgeFwd = glm::not_(edgeFwd); + edgeFwd = !edgeFwd; std::swap(outTri[0], outTri[1]); } for (const int i : {0, 1, 2, 3}) { @@ -161,7 +161,7 @@ class Partition { const vec4 nextBary = partition.vertBary[(i + 1) % 4]; for (int j = 1; j < n[i]; ++j) { partition.vertBary.push_back( - glm::mix(partition.vertBary[i], nextBary, (double)j / n[i])); + la::lerp(partition.vertBary[i], nextBary, (double)j / n[i])); } } PartitionQuad(partition.triVert, partition.vertBary, {0, 1, 2, 3}, @@ -174,7 +174,7 @@ class Partition { const vec4 nextBary = partition.vertBary[(i + 1) % 3]; for (int j = 1; j < n[i]; ++j) { partition.vertBary.push_back( - glm::mix(partition.vertBary[i], nextBary, (double)j / n[i])); + la::lerp(partition.vertBary[i], nextBary, (double)j / n[i])); } } const ivec3 edgeOffsets = {3, 3 + n[0] - 1, 3 + n[0] - 1 + n[1] - 1}; @@ -205,7 +205,7 @@ class Partition { const vec4 middleBary = partition.vertBary[edgeOffsets[0] + ns - 1]; for (int j = 1; j < nh; ++j) { partition.vertBary.push_back( - glm::mix(partition.vertBary[2], middleBary, (double)j / nh)); + la::lerp(partition.vertBary[2], middleBary, (double)j / nh)); } partition.triVert.push_back({edgeOffsets[1] - 1, 1, edgeOffsets[1]}); @@ -259,12 +259,12 @@ class Partition { // are zero, in which case a terminal triangulation is performed. static void PartitionQuad(Vec& triVert, Vec& vertBary, ivec4 cornerVerts, ivec4 edgeOffsets, - ivec4 edgeAdded, glm::bvec4 edgeFwd) { + ivec4 edgeAdded, bvec4 edgeFwd) { auto GetEdgeVert = [&](int edge, int idx) { return edgeOffsets[edge] + (edgeFwd[edge] ? 1 : -1) * idx; }; - DEBUG_ASSERT(glm::all(glm::greaterThanEqual(edgeAdded, ivec4(0))), logicErr, + DEBUG_ASSERT(la::all(la::gequal(edgeAdded, ivec4(0))), logicErr, "negative divisions!"); int corner = -1; @@ -326,7 +326,7 @@ class Partition { ivec4 newEdgeOffsets = {edgeOffsets[1], -1, GetEdgeVert(3, edgeAdded[3] + 1), edgeOffsets[0]}; ivec4 newEdgeAdded = {0, -1, 0, edgeAdded[0]}; - glm::bvec4 newEdgeFwd = {edgeFwd[1], true, edgeFwd[3], edgeFwd[0]}; + bvec4 newEdgeFwd = {edgeFwd[1], true, edgeFwd[3], edgeFwd[0]}; for (int i = 1; i < partitions; ++i) { const int cornerOffset1 = (edgeAdded[1] * i) / partitions; @@ -334,7 +334,7 @@ class Partition { edgeAdded[3] - 1 - (edgeAdded[3] * i) / partitions; const int nextOffset1 = GetEdgeVert(1, cornerOffset1 + 1); const int nextOffset3 = GetEdgeVert(3, cornerOffset3 + 1); - const int added = std::round(glm::mix( + const int added = std::round(la::lerp( (double)edgeAdded[0], (double)edgeAdded[2], (double)i / partitions)); newCornerVerts[1] = GetEdgeVert(1, cornerOffset1); @@ -346,7 +346,7 @@ class Partition { newEdgeOffsets[2] = nextOffset3; for (int j = 0; j < added; ++j) { - vertBary.push_back(glm::mix(vertBary[newCornerVerts[1]], + vertBary.push_back(la::lerp(vertBary[newCornerVerts[1]], vertBary[newCornerVerts[2]], (j + 1.0) / (added + 1.0))); } @@ -461,7 +461,7 @@ void Manifold::Impl::FillRetainedVerts(Vec& vertBary) const { for (const int i : {0, 1, 2}) { const BaryIndices indices = GetIndices(3 * tri + i); if (indices.start4 < 0) continue; // skip quad interiors - vec4 uvw(0); + vec4 uvw(0.0); uvw[indices.start4] = 1; vertBary[halfedge_[3 * tri + i].startVert] = {indices.tri, uvw}; } @@ -505,11 +505,12 @@ Vec Manifold::Impl::Subdivide( return; } const vec3 vec = vertPos_[edge.first] - vertPos_[edge.second]; - const vec4 tangent0 = - halfedgeTangent_.empty() ? vec4(0) : halfedgeTangent_[hIdx]; + const vec4 tangent0 = halfedgeTangent_.empty() + ? vec4(0.0) + : halfedgeTangent_[hIdx]; const vec4 tangent1 = halfedgeTangent_.empty() - ? vec4(0) + ? vec4(0.0) : halfedgeTangent_[halfedge_[hIdx].pairedHalfedge]; edgeAdded[i] = edgeDivisions(vec, tangent0, tangent1); }); @@ -536,7 +537,7 @@ Vec Manifold::Impl::Subdivide( int total = 0; for (int j : {0, 1, 2}) { const int added = edgeAdded[half2Edge[hIdx]]; - longest = glm::max(longest, added); + longest = la::max(longest, added); total += added; hIdx = NextHalfedge(hIdx); if (IsMarkedInsideQuad(hIdx)) { @@ -551,8 +552,7 @@ Vec Manifold::Impl::Subdivide( return extra > 0 ? (extra * (longest - thisAdded)) / longest : 0; }; - tmp[i] += - glm::max(Added(hIdx), Added(halfedge_[hIdx].pairedHalfedge)); + tmp[i] += la::max(Added(hIdx), Added(halfedge_[hIdx].pairedHalfedge)); }); edgeAdded.swap(tmp); } @@ -576,7 +576,7 @@ Vec Manifold::Impl::Subdivide( const double frac = 1.0 / (n + 1); for (int i = 0; i < n; ++i) { - vec4 uvw(0); + vec4 uvw(0.0); uvw[indices.end4] = (i + 1) * frac; uvw[indices.start4] = 1 - uvw[indices.end4]; vertBary[offset + i].uvw = uvw; @@ -625,7 +625,7 @@ Vec Manifold::Impl::Subdivide( if (halfedges[0] < 0) return; ivec4 tri3; ivec4 edgeOffsets; - glm::bvec4 edgeFwd(false); + bvec4 edgeFwd(false); for (const int i : {0, 1, 2, 3}) { if (halfedges[i] < 0) { tri3[i] = -1; @@ -675,7 +675,7 @@ Vec Manifold::Impl::Subdivide( } newVertPos[vert] = triPos * vec3(bary.uvw); } else { - mat4x3 quadPos; + mat3x4 quadPos; for (const int i : {0, 1, 2, 3}) { quadPos[i] = vertPos_[halfedge_[halfedges[i]].startVert]; } @@ -715,7 +715,7 @@ Vec Manifold::Impl::Subdivide( rel.numProp + p]; } - prop[vert * rel.numProp + p] = glm::dot(triProp, vec3(bary.uvw)); + prop[vert * rel.numProp + p] = la::dot(triProp, vec3(bary.uvw)); } else { vec4 quadProp; for (const int i : {0, 1, 2, 3}) { @@ -724,7 +724,7 @@ Vec Manifold::Impl::Subdivide( quadProp[i] = rel.properties[rel.triProperties[tri][j] * rel.numProp + p]; } - prop[vert * rel.numProp + p] = glm::dot(quadProp, bary.uvw); + prop[vert * rel.numProp + p] = la::dot(quadProp, bary.uvw); } } }); @@ -747,7 +747,7 @@ Vec Manifold::Impl::Subdivide( for (int i = 0; i < n; ++i) { for (int p = 0; p < rel.numProp; ++p) { prop[(offset + i) * rel.numProp + p] = - glm::mix(rel.properties[prop0 * rel.numProp + p], + la::lerp(rel.properties[prop0 * rel.numProp + p], rel.properties[prop1 * rel.numProp + p], (i + 1) * frac); } @@ -765,7 +765,7 @@ Vec Manifold::Impl::Subdivide( auto& rel = meshRelation_; ivec4 tri3; ivec4 edgeOffsets; - glm::bvec4 edgeFwd(true); + bvec4 edgeFwd(true); for (const int i : {0, 1, 2, 3}) { if (halfedges[i] < 0) { tri3[i] = -1; diff --git a/src/svd.h b/src/svd.h index 4bb532e47..f702b42b6 100644 --- a/src/svd.h +++ b/src/svd.h @@ -82,7 +82,7 @@ struct QR { mat3 Q, R; }; // Calculates the squared norm of the vector. -inline double Dist2(vec3 v) { return glm::dot(v, v); } +inline double Dist2(vec3 v) { return la::dot(v, v); } // For an explanation of the math see // http://pages.cs.wisc.edu/~sifakis/papers/SVD_TR1690.pdf Computing the // Singular Value Decomposition of 3 x 3 matrices with minimal branching and @@ -148,15 +148,15 @@ inline mat3 JacobiEigenAnalysis(Symmetric3x3 S) { JacobiConjugation(1, 2, 0, S, q); JacobiConjugation(2, 0, 1, S, q); } - return mat3(1.0 - 2.0 * (fma(q.y, q.y, q.z * q.z)), // - 2.0 * fma(q.x, q.y, +q.w * q.z), // - 2.0 * fma(q.x, q.z, -q.w * q.y), // - 2 * fma(q.x, q.y, -q.w * q.z), // - 1 - 2 * fma(q.x, q.x, q.z * q.z), // - 2 * fma(q.y, q.z, q.w * q.x), // - 2 * fma(q.x, q.z, q.w * q.y), // - 2 * fma(q.y, q.z, -q.w * q.x), // - 1 - 2 * fma(q.x, q.x, q.y * q.y)); + return mat3({1.0 - 2.0 * (fma(q.y, q.y, q.z * q.z)), // + 2.0 * fma(q.x, q.y, +q.w * q.z), // + 2.0 * fma(q.x, q.z, -q.w * q.y)}, // + {2 * fma(q.x, q.y, -q.w * q.z), // + 1 - 2 * fma(q.x, q.x, q.z * q.z), // + 2 * fma(q.y, q.z, q.w * q.x)}, // + {2 * fma(q.x, q.z, q.w * q.y), // + 2 * fma(q.y, q.z, -q.w * q.x), // + 1 - 2 * fma(q.x, q.x, q.y * q.y)}); } // Implementation of Algorithm 3 inline void SortSingularValues(mat3& B, mat3& V) { @@ -287,12 +287,12 @@ struct SVDSet { }; /** - * Returns the Singular Value Decomposition of A: A = U * S * glm::transpose(V). + * Returns the Singular Value Decomposition of A: A = U * S * la::transpose(V). * * @param A The matrix to decompose. */ inline SVDSet SVD(mat3 A) { - mat3 V = JacobiEigenAnalysis(glm::transpose(A) * A); + mat3 V = JacobiEigenAnalysis(la::transpose(A) * A); auto B = A * V; SortSingularValues(B, V); QR qr = QRDecomposition(B); diff --git a/src/tri_dist.h b/src/tri_dist.h index b3d7b19f0..4e83ca525 100644 --- a/src/tri_dist.h +++ b/src/tri_dist.h @@ -15,7 +15,6 @@ #pragma once #include -#include #include "manifold/common.h" @@ -43,11 +42,11 @@ inline void EdgeEdgeDist(vec3& x, vec3& y, // closest points const vec3& b) // seg 2 origin, vector { const vec3 T = q - p; - const auto ADotA = glm::dot(a, a); - const auto BDotB = glm::dot(b, b); - const auto ADotB = glm::dot(a, b); - const auto ADotT = glm::dot(a, T); - const auto BDotT = glm::dot(b, T); + const auto ADotA = la::dot(a, a); + const auto BDotB = la::dot(b, b); + const auto ADotB = la::dot(a, b); + const auto ADotT = la::dot(a, T); + const auto BDotT = la::dot(b, T); // t parameterizes ray (p, a) // u parameterizes ray (q, b) @@ -57,7 +56,7 @@ inline void EdgeEdgeDist(vec3& x, vec3& y, // closest points double t; // We will clamp result so t is on the segment (p, a) t = Denom != 0.0 - ? glm::clamp((ADotT * BDotB - BDotT * ADotB) / Denom, 0.0, 1.0) + ? la::clamp((ADotT * BDotB - BDotT * ADotB) / Denom, 0.0, 1.0) : 0.0; // find u for point on ray (q, b) closest to point at t @@ -69,14 +68,14 @@ inline void EdgeEdgeDist(vec3& x, vec3& y, // closest points // otherwise, clamp u, recompute and clamp t if (u < 0.0) { u = 0.0; - t = ADotA != 0.0 ? glm::clamp(ADotT / ADotA, 0.0, 1.0) : 0.0; + t = ADotA != 0.0 ? la::clamp(ADotT / ADotA, 0.0, 1.0) : 0.0; } else if (u > 1.0) { u = 1.0; - t = ADotA != 0.0 ? glm::clamp((ADotB + ADotT) / ADotA, 0.0, 1.0) : 0.0; + t = ADotA != 0.0 ? la::clamp((ADotB + ADotT) / ADotA, 0.0, 1.0) : 0.0; } } else { u = 0.0; - t = ADotA != 0.0 ? glm::clamp(ADotT / ADotA, 0.0, 1.0) : 0.0; + t = ADotA != 0.0 ? la::clamp(ADotT / ADotA, 0.0, 1.0) : 0.0; } x = p + a * t; y = q + b * u; @@ -115,7 +114,7 @@ inline auto DistanceTriangleTriangleSquared(const std::array& p, vec3 cq; EdgeEdgeDist(cp, cq, p[i], Sv[i], q[j], Tv[j]); const vec3 V = cq - cp; - const auto dd = glm::dot(V, V); + const auto dd = la::dot(V, V); if (dd <= mindd) { mindd = dd; @@ -123,14 +122,14 @@ inline auto DistanceTriangleTriangleSquared(const std::array& p, uint32_t id = i + 2; if (id >= 3) id -= 3; vec3 Z = p[id] - cp; - auto a = glm::dot(Z, V); + auto a = la::dot(Z, V); id = j + 2; if (id >= 3) id -= 3; Z = q[id] - cq; - auto b = glm::dot(Z, V); + auto b = la::dot(Z, V); if ((a <= 0.0) && (b >= 0.0)) { - return glm::dot(V, V); + return la::dot(V, V); }; if (a <= 0.0) @@ -143,12 +142,12 @@ inline auto DistanceTriangleTriangleSquared(const std::array& p, } } - vec3 Sn = glm::cross(Sv[0], Sv[1]); - auto Snl = glm::dot(Sn, Sn); + vec3 Sn = la::cross(Sv[0], Sv[1]); + auto Snl = la::dot(Sn, Sn); if (Snl > 1e-15) { - const vec3 Tp(glm::dot(p[0] - q[0], Sn), glm::dot(p[0] - q[1], Sn), - glm::dot(p[0] - q[2], Sn)); + const vec3 Tp(la::dot(p[0] - q[0], Sn), la::dot(p[0] - q[1], Sn), + la::dot(p[0] - q[2], Sn)); int index = -1; if ((Tp[0] > 0.0) && (Tp[1] > 0.0) && (Tp[2] > 0.0)) { @@ -165,29 +164,29 @@ inline auto DistanceTriangleTriangleSquared(const std::array& p, const vec3& qIndex = q[index]; vec3 V = qIndex - p[0]; - vec3 Z = glm::cross(Sn, Sv[0]); - if (glm::dot(V, Z) > 0.0) { + vec3 Z = la::cross(Sn, Sv[0]); + if (la::dot(V, Z) > 0.0) { V = qIndex - p[1]; - Z = glm::cross(Sn, Sv[1]); - if (glm::dot(V, Z) > 0.0) { + Z = la::cross(Sn, Sv[1]); + if (la::dot(V, Z) > 0.0) { V = qIndex - p[2]; - Z = glm::cross(Sn, Sv[2]); - if (glm::dot(V, Z) > 0.0) { + Z = la::cross(Sn, Sv[2]); + if (la::dot(V, Z) > 0.0) { vec3 cp = qIndex + Sn * Tp[index] / Snl; vec3 cq = qIndex; - return glm::dot(cp - cq, cp - cq); + return la::dot(cp - cq, cp - cq); } } } } } - vec3 Tn = glm::cross(Tv[0], Tv[1]); - auto Tnl = glm::dot(Tn, Tn); + vec3 Tn = la::cross(Tv[0], Tv[1]); + auto Tnl = la::dot(Tn, Tn); if (Tnl > 1e-15) { - const vec3 Sp(glm::dot(q[0] - p[0], Tn), glm::dot(q[0] - p[1], Tn), - glm::dot(q[0] - p[2], Tn)); + const vec3 Sp(la::dot(q[0] - p[0], Tn), la::dot(q[0] - p[1], Tn), + la::dot(q[0] - p[2], Tn)); int index = -1; if ((Sp[0] > 0.0) && (Sp[1] > 0.0) && (Sp[2] > 0.0)) { @@ -204,17 +203,17 @@ inline auto DistanceTriangleTriangleSquared(const std::array& p, const vec3& pIndex = p[index]; vec3 V = pIndex - q[0]; - vec3 Z = glm::cross(Tn, Tv[0]); - if (glm::dot(V, Z) > 0.0) { + vec3 Z = la::cross(Tn, Tv[0]); + if (la::dot(V, Z) > 0.0) { V = pIndex - q[1]; - Z = glm::cross(Tn, Tv[1]); - if (glm::dot(V, Z) > 0.0) { + Z = la::cross(Tn, Tv[1]); + if (la::dot(V, Z) > 0.0) { V = pIndex - q[2]; - Z = glm::cross(Tn, Tv[2]); - if (glm::dot(V, Z) > 0.0) { + Z = la::cross(Tn, Tv[2]); + if (la::dot(V, Z) > 0.0) { vec3 cp = pIndex; vec3 cq = pIndex + Tn * Sp[index] / Tnl; - return glm::dot(cp - cq, cp - cq); + return la::dot(cp - cq, cp - cq); } } } diff --git a/src/utils.h b/src/utils.h index f40707a75..b7e43c955 100644 --- a/src/utils.h +++ b/src/utils.h @@ -223,7 +223,7 @@ inline int CCW(vec2 p0, vec2 p1, vec2 p2, double tol) { vec2 v1 = p1 - p0; vec2 v2 = p2 - p0; double area = fma(v1.x, v2.y, -v1.y * v2.x); - double base2 = glm::max(glm::dot(v1, v1), glm::dot(v2, v2)); + double base2 = la::max(la::dot(v1, v1), la::dot(v2, v2)); if (area * area * 4 <= base2 * tol * tol) return 0; else @@ -236,14 +236,20 @@ inline int CCW(vec2 p0, vec2 p1, vec2 p2, double tol) { * * @param up The vector to be turned to point upwards. Length does not matter. */ -inline mat4x3 RotateUp(vec3 up) { - up = glm::normalize(up); - vec3 axis = glm::cross(up, {0, 0, 1}); - double angle = glm::asin(glm::length(axis)); - if (glm::dot(up, {0, 0, 1}) < 0) angle = glm::pi() - angle; - return mat4x3(glm::rotate(mat4(1), angle, axis)); +inline mat3x4 RotateUp(vec3 up) { + up = la::normalize(up); + const vec3 axis = la::cross(up, {0, 0, 1}); + double angle = la::asin(la::length(axis)); + if (la::dot(up, {0, 0, 1}) < 0) angle = kPi - angle; + const quat q = la::rotation_quat(la::normalize(axis), angle); + return mat3x4(la::qmat(q), vec3()); } +inline mat4 Mat4(mat3x4 a) { + return mat4({a[0], 0}, {a[1], 0}, {a[2], 0}, {a[3], 1}); +} +inline mat3 Mat3(mat2x3 a) { return mat3({a[0], 0}, {a[1], 0}, {a[2], 1}); } + /** @} */ /** @defgroup Debug @@ -255,34 +261,38 @@ inline mat4x3 RotateUp(vec3 up) { */ #ifdef MANIFOLD_DEBUG -template -inline std::ostream& operator<<(std::ostream& stream, const glm::tvec2& v) { - return stream << "x = " << v.x << ", y = " << v.y; +template +std::ostream& operator<<(std::ostream& out, const la::vec& v) { + return out << '{' << v[0] << '}'; } - -template -inline std::ostream& operator<<(std::ostream& stream, const glm::tvec3& v) { - return stream << "x = " << v.x << ", y = " << v.y << ", z = " << v.z; +template +std::ostream& operator<<(std::ostream& out, const la::vec& v) { + return out << '{' << v[0] << ',' << v[1] << '}'; } - -template -inline std::ostream& operator<<(std::ostream& stream, const glm::tvec4& v) { - return stream << "x = " << v.x << ", y = " << v.y << ", z = " << v.z - << ", w = " << v.w; +template +std::ostream& operator<<(std::ostream& out, const la::vec& v) { + return out << '{' << v[0] << ',' << v[1] << ',' << v[2] << '}'; } - -inline std::ostream& operator<<(std::ostream& stream, const mat3& mat) { - mat3 tam = glm::transpose(mat); - return stream << tam[0] << std::endl - << tam[1] << std::endl - << tam[2] << std::endl; +template +std::ostream& operator<<(std::ostream& out, const la::vec& v) { + return out << '{' << v[0] << ',' << v[1] << ',' << v[2] << ',' << v[3] << '}'; } -inline std::ostream& operator<<(std::ostream& stream, const mat4x3& mat) { - mat3x4 tam = glm::transpose(mat); - return stream << tam[0] << std::endl - << tam[1] << std::endl - << tam[2] << std::endl; +template +std::ostream& operator<<(std::ostream& out, const la::mat& m) { + return out << '{' << m[0] << '}'; +} +template +std::ostream& operator<<(std::ostream& out, const la::mat& m) { + return out << '{' << m[0] << ',' << m[1] << '}'; +} +template +std::ostream& operator<<(std::ostream& out, const la::mat& m) { + return out << '{' << m[0] << ',' << m[1] << ',' << m[2] << '}'; +} +template +std::ostream& operator<<(std::ostream& out, const la::mat& m) { + return out << '{' << m[0] << ',' << m[1] << ',' << m[2] << ',' << m[3] << '}'; } inline std::ostream& operator<<(std::ostream& stream, const Box& box) { diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index eb95d43ba..2a64fad7e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -45,18 +45,19 @@ if(NOT gtest_FOUND) FetchContent_MakeAvailable(googletest) endif() +# put fast/simple tests files first to run earlier set( SOURCE_FILES - boolean_complex_test.cpp - boolean_test.cpp - hull_test.cpp - manifold_test.cpp + test_main.cpp polygon_test.cpp properties_test.cpp - samples_test.cpp + manifold_test.cpp + boolean_test.cpp sdf_test.cpp smooth_test.cpp - test_main.cpp + hull_test.cpp + samples_test.cpp + boolean_complex_test.cpp ) if(MANIFOLD_CROSS_SECTION) diff --git a/test/boolean_complex_test.cpp b/test/boolean_complex_test.cpp index 7c5674cd0..5f99e88eb 100644 --- a/test/boolean_complex_test.cpp +++ b/test/boolean_complex_test.cpp @@ -142,7 +142,7 @@ TEST(BooleanComplex, Cylinders) { Manifold m1; for (auto& array : arrays1) { - mat4x3 mat; + mat3x4 mat; for (const int i : {0, 1, 2, 3}) { for (const int j : {0, 1, 2}) { mat[i][j] = array[j * 4 + i]; @@ -153,7 +153,7 @@ TEST(BooleanComplex, Cylinders) { Manifold m2; for (auto& array : arrays2) { - mat4x3 mat; + mat3x4 mat; for (const int i : {0, 1, 2, 3}) { for (const int j : {0, 1, 2}) { mat[i][j] = array[j * 4 + i]; @@ -240,9 +240,8 @@ TEST(BooleanComplex, Close) { } auto prop = result.GetProperties(); const double tol = 0.004; - EXPECT_NEAR(prop.volume, (4.0 / 3.0) * glm::pi() * r * r * r, - tol * r * r * r); - EXPECT_NEAR(prop.surfaceArea, 4 * glm::pi() * r * r, tol * r * r); + EXPECT_NEAR(prop.volume, (4.0 / 3.0) * kPi * r * r * r, tol * r * r * r); + EXPECT_NEAR(prop.surfaceArea, 4 * kPi * r * r, tol * r * r); #ifdef MANIFOLD_EXPORT if (options.exportModels) ExportMesh("close.glb", result.GetMeshGL(), {}); @@ -252,17 +251,13 @@ TEST(BooleanComplex, Close) { } TEST(BooleanComplex, BooleanVolumes) { - mat4 m = glm::translate(mat4(1.0), vec3(1.0)); - // Define solids which volumes are easy to compute w/ bit arithmetics: // m1, m2, m4 are unique, non intersecting "bits" (of volume 1, 2, 4) // m3 = m1 + m2 // m7 = m1 + m2 + m3 auto m1 = Manifold::Cube({1, 1, 1}); - auto m2 = Manifold::Cube({2, 1, 1}).Transform( - mat4x3(glm::translate(mat4(1.0), vec3(1.0, 0, 0)))); - auto m4 = Manifold::Cube({4, 1, 1}).Transform( - mat4x3(glm::translate(mat4(1.0), vec3(3.0, 0, 0)))); + auto m2 = Manifold::Cube({2, 1, 1}).Translate({1, 0, 0}); + auto m4 = Manifold::Cube({4, 1, 1}).Translate({3, 0, 0}); auto m3 = Manifold::Cube({3, 1, 1}); auto m7 = Manifold::Cube({7, 1, 1}); @@ -283,7 +278,7 @@ TEST(BooleanComplex, Spiral) { const int d = 2; std::function spiral = [&](const int rec, const double r, const double add) { - const double rot = 360.0 / (glm::pi() * r * 2) * d; + const double rot = 360.0 / (kPi * r * 2) * d; const double rNext = r + add / 360 * rot; const Manifold cube = Manifold::Cube(vec3(1), true).Translate({0, r, 0}); @@ -301,9 +296,9 @@ TEST(BooleanComplex, Sweep) { // generate the minimum equivalent positive angle auto minPosAngle = [](double angle) { - double div = angle / glm::two_pi(); + double div = angle / kTwoPi; double wholeDiv = floor(div); - return angle - wholeDiv * glm::two_pi(); + return angle - wholeDiv * kTwoPi; }; // calculate determinant @@ -318,7 +313,7 @@ TEST(BooleanComplex, Sweep) { std::vector arcPoints; for (int i = 0; i < numberOfArcPoints; i++) { - double angle = i * glm::pi() / numberOfArcPoints; + double angle = i * kPi / numberOfArcPoints; double y = arcCenterPoint.y - cos(angle) * filletRadius; double x = arcCenterPoint.x + sin(angle) * filletRadius; arcPoints.push_back(vec2(x, y)); @@ -349,8 +344,7 @@ TEST(BooleanComplex, Sweep) { totalAngle = posEndAngle - startAngle; } - int nSegments = - ceil(totalAngle / glm::two_pi() * nSegmentsPerRotation + 1); + int nSegments = ceil(totalAngle / kTwoPi * nSegmentsPerRotation + 1); if (nSegments < 2) { nSegments = 2; } @@ -389,7 +383,7 @@ TEST(BooleanComplex, Sweep) { Manifold::Extrude(profile.ToPolygons(), distance) .Rotate(90, 0, -90) .Translate(vec3(distance, 0, 0)) - .Rotate(0, 0, angle * 180 / glm::pi()) + .Rotate(0, 0, angle * 180 / kPi) .Translate(vec3(p1.x, p1.y, 0)); std::vector result; diff --git a/test/cross_section_test.cpp b/test/cross_section_test.cpp index 004433488..7ff81a6df 100644 --- a/test/cross_section_test.cpp +++ b/test/cross_section_test.cpp @@ -43,7 +43,7 @@ TEST(CrossSection, MirrorUnion) { #endif EXPECT_FLOAT_EQ(2.5 * a.Area(), cross.Area()); - EXPECT_TRUE(a.Mirror(vec2(0)).IsEmpty()); + EXPECT_TRUE(a.Mirror(vec2(0.0)).IsEmpty()); } TEST(CrossSection, RoundOffset) { @@ -88,17 +88,17 @@ TEST(CrossSection, Transform) { auto sq = CrossSection::Square({10., 10.}); auto a = sq.Rotate(45).Scale({2, 3}).Translate({4, 5}); - mat3 trans(1.0, 0.0, 0.0, // - 0.0, 1.0, 0.0, // - 4.0, 5.0, 1.0); - mat3 rot(cosd(45), sind(45), 0.0, // - -sind(45), cosd(45), 0.0, // - 0.0, 0.0, 1.0); - mat3 scale(2.0, 0.0, 0.0, // - 0.0, 3.0, 0.0, // - 0.0, 0.0, 1.0); - - auto b = sq.Transform(mat3x2(trans * scale * rot)); + mat3 trans({1.0, 0.0, 0.0}, // + {0.0, 1.0, 0.0}, // + {4.0, 5.0, 1.0}); + mat3 rot({cosd(45), sind(45), 0.0}, // + {-sind(45), cosd(45), 0.0}, // + {0.0, 0.0, 1.0}); + mat3 scale({2.0, 0.0, 0.0}, // + {0.0, 3.0, 0.0}, // + {0.0, 0.0, 1.0}); + + auto b = sq.Transform(mat2x3(trans * scale * rot)); auto b_copy = CrossSection(b); auto ex_b = Manifold::Extrude(b.ToPolygons(), 1.).GetMeshGL(); diff --git a/test/hull_test.cpp b/test/hull_test.cpp index 2bc4cc892..08106399c 100644 --- a/test/hull_test.cpp +++ b/test/hull_test.cpp @@ -39,7 +39,7 @@ bool isMeshConvex(Manifold hullManifold, double epsilon = 0.0000001) { vec3 v2 = mesh.GetVertPos(tri[2]); // Compute the normal of the triangle - vec3 normal = glm::normalize(glm::cross(v1 - v0, v2 - v0)); + vec3 normal = la::normalize(la::cross(v1 - v0, v2 - v0)); // Check all other vertices for (size_t i = 0; i < numVert; ++i) { @@ -50,7 +50,7 @@ bool isMeshConvex(Manifold hullManifold, double epsilon = 0.0000001) { vec3 v = mesh.GetVertPos(i); // Compute the signed distance from the plane - double distance = glm::dot(normal, v - v0); + double distance = la::dot(normal, v - v0); // If any vertex lies on the opposite side of the normal direction if (distance > epsilon) { diff --git a/test/manifold_test.cpp b/test/manifold_test.cpp index adb6c6c6b..fbc0c922f 100644 --- a/test/manifold_test.cpp +++ b/test/manifold_test.cpp @@ -213,8 +213,8 @@ TEST(Manifold, Revolve) { vug = Manifold::Revolve(rotatedPolys, 48); EXPECT_EQ(vug.Genus(), -1); auto prop = vug.GetProperties(); - EXPECT_NEAR(prop.volume, 14.0 * glm::pi(), 0.2); - EXPECT_NEAR(prop.surfaceArea, 30.0 * glm::pi(), 0.2); + EXPECT_NEAR(prop.volume, 14.0 * kPi, 0.2); + EXPECT_NEAR(prop.surfaceArea, 30.0 * kPi, 0.2); } } @@ -223,8 +223,8 @@ TEST(Manifold, Revolve2) { Manifold donutHole = Manifold::Revolve(polys, 48); EXPECT_EQ(donutHole.Genus(), 0); auto prop = donutHole.GetProperties(); - EXPECT_NEAR(prop.volume, 48.0 * glm::pi(), 1.0); - EXPECT_NEAR(prop.surfaceArea, 96.0 * glm::pi(), 1.0); + EXPECT_NEAR(prop.volume, 48.0 * kPi, 1.0); + EXPECT_NEAR(prop.surfaceArea, 96.0 * kPi, 1.0); } #ifdef MANIFOLD_CROSS_SECTION @@ -232,8 +232,8 @@ TEST(Manifold, Revolve3) { CrossSection circle = CrossSection::Circle(1, 32); Manifold sphere = Manifold::Revolve(circle.ToPolygons(), 32); auto prop = sphere.GetProperties(); - EXPECT_NEAR(prop.volume, 4.0 / 3.0 * glm::pi(), 0.1); - EXPECT_NEAR(prop.surfaceArea, 4 * glm::pi(), 0.15); + EXPECT_NEAR(prop.volume, 4.0 / 3.0 * kPi, 0.1); + EXPECT_NEAR(prop.surfaceArea, 4 * kPi, 0.15); } #endif @@ -247,10 +247,9 @@ TEST(Manifold, PartialRevolveOnYAxis) { revolute = Manifold::Revolve(rotatedPolys, 48, 180); EXPECT_EQ(revolute.Genus(), 1); auto prop = revolute.GetProperties(); - EXPECT_NEAR(prop.volume, 24.0 * glm::pi(), 1.0); + EXPECT_NEAR(prop.volume, 24.0 * kPi, 1.0); EXPECT_NEAR(prop.surfaceArea, - 48.0 * glm::pi() + 4.0 * 4.0 * 2.0 - 2.0 * 2.0 * 2.0, - 1.0); + 48.0 * kPi + 4.0 * 4.0 * 2.0 - 2.0 * 2.0 * 2.0, 1.0); } } @@ -291,7 +290,7 @@ TEST(Manifold, Warp2) { Manifold shape = Manifold::Extrude(circle.ToPolygons(), 2, 10).Warp([](vec3& v) { int nSegments = 10; - double angleStep = 2.0 / 3.0 * glm::pi() / nSegments; + double angleStep = 2.0 / 3.0 * kPi / nSegments; int zIndex = nSegments - 1 - std::round(v.z); double angle = zIndex * angleStep; v.z = v.y; @@ -422,20 +421,20 @@ TEST(Manifold, Transform) { Manifold cube2 = cube; cube = cube.Rotate(30, 40, 50).Scale({6, 5, 4}).Translate({1, 2, 3}); - mat3 rX(1.0, 0.0, 0.0, // - 0.0, cosd(30), sind(30), // - 0.0, -sind(30), cosd(30)); - mat3 rY(cosd(40), 0.0, -sind(40), // - 0.0, 1.0, 0.0, // - sind(40), 0.0, cosd(40)); - mat3 rZ(cosd(50), sind(50), 0.0, // - -sind(50), cosd(50), 0.0, // - 0.0, 0.0, 1.0); - mat3 s = mat3(1.0); + mat3 rX({1.0, 0.0, 0.0}, // + {0.0, cosd(30), sind(30)}, // + {0.0, -sind(30), cosd(30)}); + mat3 rY({cosd(40), 0.0, -sind(40)}, // + {0.0, 1.0, 0.0}, // + {sind(40), 0.0, cosd(40)}); + mat3 rZ({cosd(50), sind(50), 0.0}, // + {-sind(50), cosd(50), 0.0}, // + {0.0, 0.0, 1.0}); + mat3 s; s[0][0] = 6; s[1][1] = 5; s[2][2] = 4; - mat4x3 transform = mat4x3(s * rZ * rY * rX); + mat3x4 transform = mat3x4(s * rZ * rY * rX, vec3(0.0)); transform[3] = vec3(1, 2, 3); cube2 = cube2.Transform(transform); @@ -611,7 +610,7 @@ TEST(Manifold, MirrorUnion) { auto vol_a = a.GetProperties().volume; EXPECT_FLOAT_EQ(vol_a * 2.75, result.GetProperties().volume); - EXPECT_TRUE(a.Mirror(vec3(0)).IsEmpty()); + EXPECT_TRUE(a.Mirror(vec3(0.0)).IsEmpty()); } TEST(Manifold, MirrorUnion2) { diff --git a/test/manifoldc_test.cpp b/test/manifoldc_test.cpp index ef20ae809..c8e6dfda9 100644 --- a/test/manifoldc_test.cpp +++ b/test/manifoldc_test.cpp @@ -143,12 +143,13 @@ TEST(CBIND, level_set) { double a = context[0] * context[1]; double b = context[0] * context[2]; double c = context[0] * context[3]; - double s = 4.0 * glm::pi() * + constexpr double kPi = 3.14159265358979323846264338327950288; + double s = 4.0 * kPi * std::pow(((std::pow(a * b, 1.6) + std::pow(a * c, 1.6) + std::pow(b * c, 1.6)) / 3.0), 1.0 / 1.6); - double v = 4.0 * glm::pi() / 3.0 * a * b * c; + double v = 4.0 * kPi / 3.0 * a * b * c; // Numerical calculations for volume and surface area ManifoldProperties sdf_props = manifold_get_properties(sdf_man); @@ -174,7 +175,7 @@ TEST(CBIND, properties) { void *) = [](double *new_prop, ManifoldVec3 position, const double *old_prop, void *ctx) { new_prop[0] = - glm::sqrt(glm::sqrt(position.x * position.x + position.y * position.y) + + std::sqrt(std::sqrt(position.x * position.x + position.y * position.y) + position.z * position.z) * 5.0; }; @@ -184,7 +185,7 @@ TEST(CBIND, properties) { void *) = [](double *new_prop, ManifoldVec3 position, const double *old_prop, void *ctx) { new_prop[0] = - glm::sqrt(glm::sqrt(position.x * position.x + position.y * position.y) + + std::sqrt(std::sqrt(position.x * position.x + position.y * position.y) + position.z * position.z) * ((double *)ctx)[0]; }; diff --git a/test/polygon_test.cpp b/test/polygon_test.cpp index f433cc656..ba11fd293 100644 --- a/test/polygon_test.cpp +++ b/test/polygon_test.cpp @@ -27,7 +27,7 @@ using namespace manifold; Polygons Turn180(Polygons polys) { for (SimplePolygon &poly : polys) { for (vec2 &vert : poly) { - vert *= -1; + vert *= -1.0; } } return polys; diff --git a/test/samples_test.cpp b/test/samples_test.cpp index 170957da4..dc1dd6548 100644 --- a/test/samples_test.cpp +++ b/test/samples_test.cpp @@ -114,7 +114,7 @@ TEST(Samples, Scallop) { const vec3 red(1, 0, 0); const vec3 blue(0, 0, 1); const double limit = 15; - vec3 color = glm::mix(blue, red, glm::smoothstep(-limit, limit, curvature)); + vec3 color = la::lerp(blue, red, smoothstep(-limit, limit, curvature)); for (const int i : {0, 1, 2}) { newProp[i] = color[i]; } diff --git a/test/sdf_test.cpp b/test/sdf_test.cpp index bdb54b622..f5fe4a1b7 100644 --- a/test/sdf_test.cpp +++ b/test/sdf_test.cpp @@ -29,7 +29,7 @@ struct CubeVoid { struct Layers { double operator()(vec3 p) const { - int a = glm::mod(std::round(2 * p.z), 4.0); + int a = std::fmod(std::round(2 * p.z), 4.0); return a == 0 ? 1 : (a == 2 ? -1 : 0); } }; @@ -37,8 +37,8 @@ struct Layers { TEST(SDF, SphereShell) { Manifold sphere = Manifold::LevelSet( [](vec3 pos) { - const double r = glm::length(pos); - return glm::min(1 - r, r - 0.995f); + const double r = la::length(pos); + return la::min(1 - r, r - 0.995f); }, {vec3(-1.1), vec3(1.1)}, 0.01, 0, 0.0001); @@ -141,7 +141,7 @@ TEST(SDF, Surface) { TEST(SDF, Resize) { const double size = 20; - Manifold layers = Manifold::LevelSet(Layers(), {vec3(0), vec3(size)}, 1); + Manifold layers = Manifold::LevelSet(Layers(), {vec3(0.0), vec3(size)}, 1); #ifdef MANIFOLD_EXPORT if (options.exportModels) ExportMesh("layers.gltf", layers.GetMeshGL(), {}); #endif @@ -153,10 +153,10 @@ TEST(SDF, Resize) { TEST(SDF, SineSurface) { Manifold surface = Manifold::LevelSet( [](vec3 p) { - double mid = glm::sin(p.x) + glm::sin(p.y); + double mid = la::sin(p.x) + la::sin(p.y); return (p.z > mid - 0.5 && p.z < mid + 0.5) ? 1.0f : -1.0f; }, - {vec3(-1.75 * glm::pi()), vec3(1.75 * glm::pi())}, 1); + {vec3(-1.75 * kPi), vec3(1.75 * kPi)}, 1); Manifold smoothed = surface.SmoothOut(180).RefineToLength(0.05); EXPECT_EQ(smoothed.Status(), Manifold::Error::NoError); @@ -170,23 +170,23 @@ TEST(SDF, SineSurface) { TEST(SDF, Blobs) { const double blend = 1; - std::vector balls = {{0, 0, 0, 2}, // - {1, 2, 3, 2}, // - {-2, 2, -2, 1}, // - {-2, -3, -2, 2}, // - {-3, -1, -3, 1}, // - {2, -3, -2, 2}, // - {-2, 3, 2, 2}, // - {-2, -3, 2, 2}, // - {1, -1, 1, -2}, // - {-4, -3, -2, 1}}; + std::vector balls = {{0, 0, 0, 2}, // + {1, 2, 3, 2}, // + {-2, 2, -2, 1}, // + {-2, -3, -2, 2}, // + {-3, -1, -3, 1}, // + {2, -3, -2, 2}, // + {-2, 3, 2, 2}, // + {-2, -3, 2, 2}, // + {1, -1, 1, -2}, // + {-4, -3, -2, 1}}; Manifold blobs = Manifold::LevelSet( [&balls, blend](vec3 p) { double d = 0; for (const auto& ball : balls) { - d += glm::sign(ball.w) * - glm::smoothstep(-blend, blend, - std::abs(ball.w) - glm::length(vec3(ball) - p)); + d += (ball.w > 0 ? 1 : -1) * + smoothstep(-blend, blend, + std::abs(ball.w) - la::length(vec3(ball) - p)); } return d; }, diff --git a/test/smooth_test.cpp b/test/smooth_test.cpp index 15a8a093d..efbdfc153 100644 --- a/test/smooth_test.cpp +++ b/test/smooth_test.cpp @@ -53,8 +53,8 @@ TEST(Smooth, RefineQuads) { .RefineToLength(0.05); EXPECT_EQ(cylinder.NumTri(), 16892); auto prop = cylinder.GetProperties(); - EXPECT_NEAR(prop.volume, 2 * glm::pi(), 0.003); - EXPECT_NEAR(prop.surfaceArea, 6 * glm::pi(), 0.004); + EXPECT_NEAR(prop.volume, 2 * kPi, 0.003); + EXPECT_NEAR(prop.surfaceArea, 6 * kPi, 0.004); const MeshGL out = cylinder.GetMeshGL(); CheckGL(out); @@ -175,7 +175,7 @@ TEST(Smooth, Precision) { // Ignore end caps. const double r2 = (std::abs(a.z) < 0.001 || std::abs(a.z - height) < 0.001) ? radius * radius - : glm::dot(a1, a1); + : la::dot(a1, a1); maxR2 = std::max(maxR2, r2); minR2 = std::min(minR2, r2); } @@ -227,8 +227,7 @@ TEST(Smooth, Manual) { 3, [](double* newProp, vec3 pos, const double* oldProp) { const vec3 red(1, 0, 0); const vec3 purple(1, 0, 1); - vec3 color = - glm::mix(purple, red, glm::smoothstep(0.0, 2.0, oldProp[0])); + vec3 color = la::lerp(purple, red, smoothstep(0.0, 2.0, oldProp[0])); for (const int i : {0, 1, 2}) newProp[i] = color[i]; }); const MeshGL out = interp.GetMeshGL(); @@ -280,7 +279,7 @@ TEST(Smooth, Csaszar) { const vec3& uvw = {0.5, 0.5, 0.0}; const double alpha = std::min(uvw[0], std::min(uvw[1], uvw[2])); options.mat.vertColor[out.triVerts[3 * tri + i]] = - glm::mix(yellow, blue, glm::smoothstep(0.0, 0.2, alpha)); + la::lerp(yellow, blue, smoothstep(0.0, 0.2, alpha)); } } ExportMesh("smoothCsaszar.glb", out, options); @@ -289,16 +288,16 @@ TEST(Smooth, Csaszar) { } vec4 CircularTangent(const vec3& tangent, const vec3& edgeVec) { - const vec3 dir = glm::normalize(tangent); + const vec3 dir = la::normalize(tangent); - double weight = std::abs(glm::dot(dir, glm::normalize(edgeVec))); + double weight = std::abs(la::dot(dir, la::normalize(edgeVec))); if (weight == 0) { weight = 1; } // Quadratic weighted bezier for circular interpolation - const vec4 bz2 = weight * vec4(dir * glm::length(edgeVec) / (2 * weight), 1); + const vec4 bz2 = weight * vec4(dir * la::length(edgeVec) / (2 * weight), 1); // Equivalent cubic weighted bezier - const vec4 bz3 = glm::mix(vec4(0, 0, 0, 1), bz2, 2 / 3.0); + const vec4 bz3 = la::lerp(vec4(0, 0, 0, 1), bz2, 2 / 3.0); // Convert from homogeneous form to geometric form return vec4(vec3(bz3) / bz3.w, bz3.w); } @@ -324,16 +323,16 @@ TEST(Smooth, Torus) { const vec3 edge = v1 - v; if (edge.z == 0) { vec3 tan(v.y, -v.x, 0); - tan *= glm::sign(glm::dot(tan, edge)); + tan *= la::dot(tan, edge) < 0 ? -1.0 : 1.0; tangent = CircularTangent(tan, edge); - } else if (std::abs(glm::determinant(mat2(vec2(v), vec2(edge)))) < + } else if (std::abs(la::determinant(mat2(vec2(v), vec2(edge)))) < kTolerance) { const double theta = std::asin(v.z); vec2 xy(v); - const double r = glm::length(xy); + const double r = la::length(xy); xy = xy / r * v.z * (r > 2 ? -1.0 : 1.0); vec3 tan(xy.x, xy.y, std::cos(theta)); - tan *= glm::sign(glm::dot(tan, edge)); + tan *= la::dot(tan, edge) < 0 ? -1.0 : 1.0; tangent = CircularTangent(tan, edge); } else { tangent = {0, 0, 0, -1}; @@ -354,8 +353,8 @@ TEST(Smooth, Torus) { vec3 v(out.vertProperties[i], out.vertProperties[i + 1], out.vertProperties[i + 2]); vec3 p(v.x, v.y, 0); - p = glm::normalize(p) * 2.0; - double r = glm::length(v - p); + p = la::normalize(p) * 2.0; + double r = la::length(v - p); ASSERT_NEAR(r, 1, 0.006); maxMeanCurvature = std::max(maxMeanCurvature, std::abs(out.vertProperties[i + 3])); @@ -376,12 +375,10 @@ TEST(Smooth, SineSurface) { Manifold surface = Manifold::LevelSet( [](vec3 p) { - double mid = glm::sin(p.x) + glm::sin(p.y); + double mid = la::sin(p.x) + la::sin(p.y); return (p.z > mid - 0.5 && p.z < mid + 0.5) ? 1.0 : -1.0; }, - {vec3(-2 * glm::pi() + 0.2), - vec3(0 * glm::pi() - 0.2)}, - 1) + {vec3(-2 * kPi + 0.2), vec3(0 * kPi - 0.2)}, 1) .AsOriginal(); Manifold smoothed = @@ -434,13 +431,13 @@ TEST(Smooth, SDF) { auto sphericalGyroid = [r](vec3 p) { const double gyroid = cos(p.x) * sin(p.y) + cos(p.y) * sin(p.z) + cos(p.z) * sin(p.x); - const double d = glm::min(0.0, r - glm::length(p)); + const double d = la::min(0.0, r - la::length(p)); return gyroid - d * d / 2; }; auto gradient = [r](vec3 pos) { - const double rad = glm::length(pos); - const double d = glm::min(0.0, r - rad) / (rad > 0 ? rad : 1); + const double rad = la::length(pos); + const double d = la::min(0.0, r - rad) / (rad > 0 ? rad : 1); const vec3 sphereGrad = d * pos; const vec3 gyroidGrad(cos(pos.z) * cos(pos.x) - sin(pos.x) * sin(pos.y), cos(pos.x) * cos(pos.y) - sin(pos.y) * sin(pos.z), @@ -465,7 +462,7 @@ TEST(Smooth, SDF) { .SetProperties( 3, [gradient](double* newProp, vec3 pos, const double* oldProp) { - const vec3 normal = -glm::normalize(gradient(pos)); + const vec3 normal = -la::normalize(gradient(pos)); for (const int i : {0, 1, 2}) newProp[i] = normal[i]; }) .SmoothByNormals(0) diff --git a/test/test.h b/test/test.h index b6c1f64bf..9c5f70ead 100644 --- a/test/test.h +++ b/test/test.h @@ -23,13 +23,6 @@ #include "manifold/meshIO.h" #endif -// somehow gcc11 + gtest 1.11.0 is unable to print ivec3 -namespace glm { -inline void PrintTo(const ivec3& point, std::ostream* os) { - *os << "(" << point.x << "," << point.y << "," << point.x << ")"; -} -} // namespace glm - using namespace manifold; struct Options { diff --git a/test/test_main.cpp b/test/test_main.cpp index cd8a4a37d..a8fc43490 100644 --- a/test/test_main.cpp +++ b/test/test_main.cpp @@ -131,7 +131,7 @@ MeshGL Csaszar() { struct GyroidSDF { double operator()(vec3 p) const { const vec3 min = p; - const vec3 max = vec3(glm::two_pi()) - p; + const vec3 max = vec3(kTwoPi) - p; const double min3 = std::min(min.x, std::min(min.y, min.z)); const double max3 = std::min(max.x, std::min(max.y, max.z)); const double bound = std::min(min3, max3); @@ -142,8 +142,8 @@ struct GyroidSDF { }; Manifold Gyroid() { - const double period = glm::two_pi(); - return Manifold::LevelSet(GyroidSDF(), {vec3(0), vec3(period)}, 0.5); + const double period = kTwoPi; + return Manifold::LevelSet(GyroidSDF(), {vec3(0.0), vec3(period)}, 0.5); } MeshGL TetGL() { @@ -181,8 +181,8 @@ MeshGL CubeSTL() { } } - const vec3 normal = glm::normalize( - glm::cross(triPos[1] - triPos[0], triPos[2] - triPos[0])); + const vec3 normal = + la::normalize(la::cross(triPos[1] - triPos[0], triPos[2] - triPos[0])); for (const int i : {0, 1, 2}) { for (const int j : {0, 1, 2}) { cube.vertProperties.push_back(triPos[i][j]); @@ -292,7 +292,7 @@ void Identical(const MeshGL& mesh1, const MeshGL& mesh2) { ASSERT_EQ(mesh1.vertProperties.size() / mesh1.numProp, mesh2.vertProperties.size() / mesh2.numProp); for (size_t i = 0; i < mesh1.vertProperties.size() / mesh1.numProp; ++i) - ASSERT_LE(glm::length(mesh1.GetVertPos(i) - mesh2.GetVertPos(i)), 0.0001); + ASSERT_LE(la::length(mesh1.GetVertPos(i) - mesh2.GetVertPos(i)), 0.0001); ASSERT_EQ(mesh1.triVerts.size(), mesh2.triVerts.size()); @@ -323,10 +323,11 @@ void RelatedGL(const Manifold& out, const std::vector& originals, MeshGL output = out.GetMeshGL(normalIdx); for (size_t run = 0; run < output.runOriginalID.size(); ++run) { const float* m = output.runTransform.data() + 12 * run; - const mat4x3 transform = output.runTransform.empty() - ? mat4x3(1.0) - : mat4x3(m[0], m[1], m[2], m[3], m[4], m[5], - m[6], m[7], m[8], m[9], m[10], m[11]); + const mat3x4 transform = + output.runTransform.empty() + ? Identity3x4() + : mat3x4({m[0], m[1], m[2]}, {m[3], m[4], m[5]}, {m[6], m[7], m[8]}, + {m[9], m[10], m[11]}); size_t i = 0; for (; i < originals.size(); ++i) { ASSERT_EQ(originals[i].runOriginalID.size(), 1); @@ -341,10 +342,10 @@ void RelatedGL(const Manifold& out, const std::vector& originals, } const int inTri = output.faceID.empty() ? tri : output.faceID[tri]; ASSERT_LT(inTri, inMesh.triVerts.size() / 3); - ivec3 inTriangle = {inMesh.triVerts[3 * inTri], - inMesh.triVerts[3 * inTri + 1], - inMesh.triVerts[3 * inTri + 2]}; - inTriangle *= inMesh.numProp; + ivec3 inTriangle(inMesh.triVerts[3 * inTri], + inMesh.triVerts[3 * inTri + 1], + inMesh.triVerts[3 * inTri + 2]); + inTriangle *= static_cast(inMesh.numProp); mat3 inTriPos; mat3 outTriPos; @@ -359,10 +360,10 @@ void RelatedGL(const Manifold& out, const std::vector& originals, inTriPos[j] = transform * pos; } vec3 outNormal = - glm::cross(outTriPos[1] - outTriPos[0], outTriPos[2] - outTriPos[0]); + la::cross(outTriPos[1] - outTriPos[0], outTriPos[2] - outTriPos[0]); vec3 inNormal = - glm::cross(inTriPos[1] - inTriPos[0], inTriPos[2] - inTriPos[0]); - const double area = glm::length(inNormal); + la::cross(inTriPos[1] - inTriPos[0], inTriPos[2] - inTriPos[0]); + const double area = la::length(inNormal); if (area == 0) continue; inNormal /= area; @@ -370,16 +371,15 @@ void RelatedGL(const Manifold& out, const std::vector& originals, const int vert = output.triVerts[3 * tri + j]; vec3 edges[3]; for (int k : {0, 1, 2}) edges[k] = inTriPos[k] - outTriPos[j]; - const double volume = - glm::dot(edges[0], glm::cross(edges[1], edges[2])); + const double volume = la::dot(edges[0], la::cross(edges[1], edges[2])); ASSERT_LE(volume, area * output.precision); if (checkNormals) { vec3 normal; for (int k : {0, 1, 2}) normal[k] = output.vertProperties[vert * output.numProp + 3 + k]; - ASSERT_NEAR(glm::length(normal), 1, 0.0001); - ASSERT_GT(glm::dot(normal, outNormal), 0); + ASSERT_NEAR(la::length(normal), 1, 0.0001); + ASSERT_GT(la::dot(normal, outNormal), 0); } else { for (size_t p = 3; p < inMesh.numProp; ++p) { const double propOut = @@ -393,7 +393,7 @@ void RelatedGL(const Manifold& out, const std::vector& originals, edgesP[k] = edges[k] + inNormal * inProp[k] - inNormal * propOut; } const double volumeP = - glm::dot(edgesP[0], glm::cross(edgesP[1], edgesP[2])); + la::dot(edgesP[0], la::cross(edgesP[1], edgesP[2])); ASSERT_LE(volumeP, area * output.precision); }