From fa6d6d7ad42c3f8d5ccb5edb45a429efdad49b6c Mon Sep 17 00:00:00 2001 From: K Date: Wed, 13 Dec 2023 08:46:07 -0500 Subject: [PATCH] Add WarpBatch - array process warp with single callback --- bindings/python/examples/torus_knot.py | 88 +++++++++++++++++ bindings/python/manifold3d.cpp | 88 +++++++++++++---- src/cross_section/include/cross_section.h | 3 + src/cross_section/src/cross_section.cpp | 45 ++++++--- src/manifold/include/manifold.h | 2 + src/manifold/src/impl.cpp | 9 +- src/manifold/src/impl.h | 1 + src/manifold/src/manifold.cpp | 14 +++ src/utilities/include/vec.h | 95 +----------------- src/utilities/include/vec_view.h | 112 ++++++++++++++++++++++ test/manifold_test.cpp | 17 ++++ 11 files changed, 349 insertions(+), 125 deletions(-) create mode 100644 bindings/python/examples/torus_knot.py create mode 100644 src/utilities/include/vec_view.h diff --git a/bindings/python/examples/torus_knot.py b/bindings/python/examples/torus_knot.py new file mode 100644 index 000000000..34d69bf4d --- /dev/null +++ b/bindings/python/examples/torus_knot.py @@ -0,0 +1,88 @@ +# %% +import numpy as np +from manifold3d import * + +# Creates a classic torus knot, defined as a string wrapping periodically +# around the surface of an imaginary donut. If p and q have a common +# factor then you will get multiple separate, interwoven knots. This is +# an example of using the warp() method, thus avoiding any direct +# handling of triangles. + + +def run(): + # The number of times the thread passes through the donut hole. + p = 1 + # The number of times the thread circles the donut. + q = 3 + # Radius of the interior of the imaginary donut. + majorRadius = 25 + # Radius of the small cross-section of the imaginary donut. + minorRadius = 10 + # Radius of the small cross-section of the actual object. + threadRadius = 3.75 + # Number of linear segments making up the threadRadius circle. Default is + # getCircularSegments(threadRadius). + circularSegments = -1 + # Number of segments along the length of the knot. Default makes roughly + # square facets. + linearSegments = -1 + + # These default values recreate Matlab Knot by Emmett Lalish: + # https://www.thingiverse.com/thing:7080 + + kLoops = np.gcd(p, q) + pk = p / kLoops + qk = q / kLoops + n = ( + circularSegments + if circularSegments > 2 + else get_circular_segments(threadRadius) + ) + m = linearSegments if linearSegments > 2 else n * qk * majorRadius / threadRadius + + offset = 2 + circle = CrossSection.circle(1, n).translate([offset, 0]) + + def ax_rotate(x, theta): + a, b = (x + 1) % 3, (x + 2) % 3 + s, c = np.sin(theta), np.cos(theta) + m = np.zeros((len(theta), 4, 4), dtype=np.float32) + m[:, a, a], m[:, a, b] = c, s + m[:, b, a], m[:, b, b] = -s, c + m[:, x, x], m[:, 3, 3] = 1, 1 + return m + + def func(pts): + npts = pts.shape[0] + x, y, z = pts[:, 0], pts[:, 1], pts[:, 2] + psi = qk * np.arctan2(x, y) + theta = psi * pk / qk + x1 = np.sqrt(x * x + y * y) + phi = np.arctan2(x1 - offset, z) + + v = np.zeros((npts, 4), dtype=np.float32) + v[:, 0] = threadRadius * np.cos(phi) + v[:, 2] = threadRadius * np.sin(phi) + v[:, 3] = 1 + r = majorRadius + minorRadius * np.cos(theta) + + m1 = ax_rotate(0, -np.arctan2(pk * minorRadius, qk * r)) + m1[:, 3, 0] = minorRadius + m2 = ax_rotate(1, theta) + m2[:, 3, 0] = majorRadius + m3 = ax_rotate(2, psi) + + v = v[:, None, :] @ (m1 @ m2 @ m3) + return v[:, 0, :3] + + def func_single(v): + pts = np.array(v)[None, :] + pts = func(pts) + return tuple(pts[0]) + + return circle.revolve(int(m)).warp_batch(func) + # return circle.revolve(int(m)).warp(func_single) + + +if __name__ == "__main__": + run() diff --git a/bindings/python/manifold3d.cpp b/bindings/python/manifold3d.cpp index 00f7d89ee..9c2fd7daf 100644 --- a/bindings/python/manifold3d.cpp +++ b/bindings/python/manifold3d.cpp @@ -70,31 +70,49 @@ namespace nb = nanobind; // helper to convert std::vector to numpy template nb::ndarray> to_numpy( - std::vector> const &vecvec) { + glm::vec const *begin, + glm::vec const *end) { // transfer ownership to PyObject - T *buffer = new T[vecvec.size() * N]; + size_t nvert = end - begin; + T *buffer = new T[nvert * N]; nb::capsule mem_mgr(buffer, [](void *p) noexcept { delete[](T *) p; }); - for (int i = 0; i < vecvec.size(); i++) { + for (int i = 0; i < nvert; i++) { for (int j = 0; j < N; j++) { - buffer[i * N + j] = vecvec[i][j]; + buffer[i * N + j] = begin[i][j]; } } - return {buffer, {vecvec.size(), N}, mem_mgr}; + return {buffer, {nvert, N}, mem_mgr}; +} + +// helper to convert std::vector to numpy +template +nb::ndarray> to_numpy( + std::vector> const &vec) { + return to_numpy(vec.data(), vec.data() + vec.size()); } // helper to convert numpy to std::vector template -std::vector> to_glm_vector( - nb::ndarray> const &arr) { - std::vector> out; - out.reserve(arr.shape(0)); +void to_glm_range(nb::ndarray> const &arr, + glm::vec *begin, + glm::vec *end) { + if (arr.shape(0) != end - begin) { + throw std::runtime_error( + "received numpy.shape[0]: " + std::to_string(arr.shape(0)) + + " expected: " + std::to_string(int(end - begin))); + } for (int i = 0; i < arr.shape(0); i++) { - out.emplace_back(); for (int j = 0; j < N; j++) { - out.back()[j] = arr(i, j); + begin[i][j] = arr(i, j); } } - return out; +} +// helper to convert numpy to std::vector +template +void to_glm_vector(nb::ndarray> const &arr, + std::vector> &out) { + out.resize(arr.shape(0)); + to_glm_range(arr, out.data(), out.data() + out.size()); } using namespace manifold; @@ -102,6 +120,10 @@ using namespace manifold; typedef std::tuple Float2; typedef std::tuple Float3; +using NumpyFloatNx2 = nb::ndarray>; +using NumpyFloatNx3 = nb::ndarray>; +using NumpyUintNx3 = nb::ndarray>; + template std::vector toVector(const T *arr, size_t size) { return std::vector(arr, arr + size); @@ -150,12 +172,10 @@ NB_MODULE(manifold3d, m) { "triangulate", [](std::vector>> polys) { - std::vector> polys_vec; - polys_vec.reserve(polys.size()); - for (auto &numpy : polys) { - polys_vec.push_back(to_glm_vector(numpy)); + std::vector> polys_vec(polys.size()); + for (int i = 0; i < polys.size(); i++) { + to_glm_vector(polys[i], polys_vec[i]); } - return to_numpy(Triangulate(polys_vec)); }, "Given a list polygons (each polygon shape=(N,2) dtype=float), " @@ -309,7 +329,23 @@ NB_MODULE(manifold3d, m) { "one which overlaps, but that is not checked here, so it is up to " "the user to choose their function with discretion." "\n\n" - ":param warpFunc: A function that modifies a given vertex position.") + ":param f: A function that modifies a given vertex position.") + .def( + "warp_batch", + [](Manifold &self, + const std::function &f) { + return self.WarpBatch([&f](VecView vecs) { + NumpyFloatNx3 arr = f(to_numpy(vecs.begin(), vecs.end())); + to_glm_range(arr, vecs.begin(), vecs.end()); + }); + }, + nb::arg("f"), + "Same as Manifold.warp but calls `f` with a " + "ndarray(shape=(N,3), dtype=float) and expects an ndarray " + "of the same shape and type in return. The input array can be " + "modified and returned if desired. " + "\n\n" + ":param f: A function that modifies multiple vertex positions.") .def( "set_properties", [](Manifold &self, int newNumProp, @@ -979,6 +1015,22 @@ NB_MODULE(manifold3d, m) { "intersections are not included in the result." "\n\n" ":param warpFunc: A function that modifies a given vertex position.") + .def( + "warp_batch", + [](CrossSection &self, + const std::function &f) { + return self.WarpBatch([&f](VecView vecs) { + NumpyFloatNx2 arr = f(to_numpy(vecs.begin(), vecs.end())); + to_glm_range(arr, vecs.begin(), vecs.end()); + }); + }, + nb::arg("f"), + "Same as CrossSection.warp but calls `f` with a " + "ndarray(shape=(N,2), dtype=float) and expects an ndarray " + "of the same shape and type in return. The input array can be " + "modified and returned if desired. " + "\n\n" + ":param f: A function that modifies multiple vertex positions.") .def("simplify", &CrossSection::Simplify, nb::arg("epsilon") = 1e-6, "Remove vertices from the contours in this CrossSection that are " "less than the specified distance epsilon from an imaginary line " diff --git a/src/cross_section/include/cross_section.h b/src/cross_section/include/cross_section.h index abe32bce5..9e15bb529 100644 --- a/src/cross_section/include/cross_section.h +++ b/src/cross_section/include/cross_section.h @@ -21,6 +21,7 @@ #include "glm/ext/matrix_float3x2.hpp" #include "glm/ext/vector_float2.hpp" #include "public.h" +#include "vec_view.h" namespace manifold { @@ -97,6 +98,8 @@ class CrossSection { CrossSection Mirror(const glm::vec2 ax) const; CrossSection Transform(const glm::mat3x2& m) const; CrossSection Warp(std::function warpFunc) const; + CrossSection WarpBatch( + std::function)> warpFunc) const; CrossSection Simplify(double epsilon = 1e-6) const; // Adapted from Clipper2 docs: diff --git a/src/cross_section/src/cross_section.cpp b/src/cross_section/src/cross_section.cpp index 9e3065942..658587255 100644 --- a/src/cross_section/src/cross_section.cpp +++ b/src/cross_section/src/cross_section.cpp @@ -567,21 +567,42 @@ CrossSection CrossSection::Transform(const glm::mat3x2& m) const { */ CrossSection CrossSection::Warp( std::function warpFunc) const { - auto paths = GetPaths(); - auto warped = C2::PathsD(); - warped.reserve(paths->paths_.size()); - for (auto path : paths->paths_) { - auto sz = path.size(); - auto s = C2::PathD(sz); - for (int i = 0; i < sz; ++i) { - auto v = v2_of_pd(path[i]); - warpFunc(v); - s[i] = v2_to_pd(v); + return WarpBatch([&warpFunc](VecView vecs) { + for (glm::vec2& p : vecs) { + warpFunc(p); + } + }); +} + +/** + * Same as CrossSection::Warp but calls warpFunc with + * a VecView which is roughly equivalent to std::span + * pointing to all vec2 elements to be modified in-place + * + * @param warpFunc A function that modifies multiple vertex positions. + */ +CrossSection CrossSection::WarpBatch( + std::function)> warpFunc) const { + std::vector tmp_verts; + C2::PathsD paths = GetPaths()->paths_; // deep copy + for (C2::PathD const& path : paths) { + for (C2::PointD const& p : path) { + tmp_verts.push_back(v2_of_pd(p)); + } + } + + warpFunc(VecView(tmp_verts.data(), tmp_verts.size())); + + auto cursor = tmp_verts.begin(); + for (C2::PathD& path : paths) { + for (C2::PointD& p : path) { + p = v2_to_pd(*cursor); + ++cursor; } - warped.push_back(s); } + return CrossSection( - shared_paths(C2::Union(warped, C2::FillRule::Positive, precision_))); + shared_paths(C2::Union(paths, C2::FillRule::Positive, precision_))); } /** diff --git a/src/manifold/include/manifold.h b/src/manifold/include/manifold.h index 283b9346e..11366d5c8 100644 --- a/src/manifold/include/manifold.h +++ b/src/manifold/include/manifold.h @@ -18,6 +18,7 @@ #include "cross_section.h" #include "public.h" +#include "vec_view.h" namespace manifold { @@ -206,6 +207,7 @@ class Manifold { Manifold Transform(const glm::mat4x3&) const; Manifold Mirror(glm::vec3) const; Manifold Warp(std::function) const; + Manifold WarpBatch(std::function)>) const; Manifold SetProperties( int, std::function) const; Manifold CalculateCurvature(int gaussianIdx, int meanIdx) const; diff --git a/src/manifold/src/impl.cpp b/src/manifold/src/impl.cpp index 5b71f4999..bf405ea9b 100644 --- a/src/manifold/src/impl.cpp +++ b/src/manifold/src/impl.cpp @@ -719,7 +719,14 @@ void Manifold::Impl::MarkFailure(Error status) { } void Manifold::Impl::Warp(std::function warpFunc) { - thrust::for_each_n(thrust::host, vertPos_.begin(), NumVert(), warpFunc); + WarpBatch([&warpFunc](VecView vecs) { + thrust::for_each(thrust::host, vecs.begin(), vecs.end(), warpFunc); + }); +} + +void Manifold::Impl::WarpBatch( + std::function)> warpFunc) { + warpFunc(vertPos_.view()); CalculateBBox(); if (!IsFinite()) { MarkFailure(Error::NonFiniteVertex); diff --git a/src/manifold/src/impl.h b/src/manifold/src/impl.h index 4d12b637c..36e37fcc0 100644 --- a/src/manifold/src/impl.h +++ b/src/manifold/src/impl.h @@ -76,6 +76,7 @@ struct Manifold::Impl { void Update(); void MarkFailure(Error status); void Warp(std::function warpFunc); + void WarpBatch(std::function)> warpFunc); Impl Transform(const glm::mat4x3& transform) const; SparseIndices EdgeCollisions(const Impl& B, bool inverted = false) const; SparseIndices VertexCollisionsZ(VecView vertsIn, diff --git a/src/manifold/src/manifold.cpp b/src/manifold/src/manifold.cpp index 124a67f86..841339aac 100644 --- a/src/manifold/src/manifold.cpp +++ b/src/manifold/src/manifold.cpp @@ -569,6 +569,20 @@ Manifold Manifold::Warp(std::function warpFunc) const { return Manifold(std::make_shared(pImpl)); } +/** + * Same as Manifold::Warp but calls warpFunc with with + * a VecView which is roughly equivalent to std::span + * pointing to all vec3 elements to be modified in-place + * + * @param warpFunc A function that modifies multiple vertex positions. + */ +Manifold Manifold::WarpBatch( + std::function)> warpFunc) const { + auto pImpl = std::make_shared(*GetCsgLeafNode().GetImpl()); + pImpl->WarpBatch(warpFunc); + return Manifold(std::make_shared(pImpl)); +} + /** * Create a new copy of this manifold with updated vertex properties by * supplying a function that takes the existing position and properties as diff --git a/src/utilities/include/vec.h b/src/utilities/include/vec.h index de84724f2..f5170b054 100644 --- a/src/utilities/include/vec.h +++ b/src/utilities/include/vec.h @@ -24,6 +24,7 @@ // #include "optional_assert.h" #include "par.h" #include "public.h" +#include "vec_view.h" namespace manifold { @@ -33,100 +34,6 @@ namespace manifold { template class Vec; -/** - * View for Vec, can perform offset operation. - * This will be invalidated when the original vector is dropped or changes - * length. - */ -template -class VecView { - public: - using Iter = T *; - using IterC = const T *; - - VecView(T *ptr_, int size_) : ptr_(ptr_), size_(size_) {} - - VecView(const VecView &other) { - ptr_ = other.ptr_; - size_ = other.size_; - } - - VecView &operator=(const VecView &other) { - ptr_ = other.ptr_; - size_ = other.size_; - return *this; - } - - // allows conversion to a const VecView - operator VecView() const { return {ptr_, size_}; } - - inline const T &operator[](int i) const { - if (i < 0 || i >= size_) throw std::out_of_range("Vec out of range"); - return ptr_[i]; - } - - inline T &operator[](int i) { - if (i < 0 || i >= size_) throw std::out_of_range("Vec out of range"); - return ptr_[i]; - } - - IterC cbegin() const { return ptr_; } - IterC cend() const { return ptr_ + size_; } - - IterC begin() const { return cbegin(); } - IterC end() const { return cend(); } - - Iter begin() { return ptr_; } - Iter end() { return ptr_ + size_; } - - const T &front() const { - if (size_ == 0) - throw std::out_of_range("attempt to take the front of an empty vector"); - return ptr_[0]; - } - - const T &back() const { - if (size_ == 0) - throw std::out_of_range("attempt to take the back of an empty vector"); - return ptr_[size_ - 1]; - } - - T &front() { - if (size_ == 0) - throw std::out_of_range("attempt to take the front of an empty vector"); - return ptr_[0]; - } - - T &back() { - if (size_ == 0) - throw std::out_of_range("attempt to take the back of an empty vector"); - return ptr_[size_ - 1]; - } - - int size() const { return size_; } - - bool empty() const { return size_ == 0; } - -#ifdef MANIFOLD_DEBUG - void Dump() { - std::cout << "Vec = " << std::endl; - for (int i = 0; i < size(); ++i) { - std::cout << i << ", " << ptr_[i] << ", " << std::endl; - } - std::cout << std::endl; - } -#endif - - protected: - T *ptr_ = nullptr; - int size_ = 0; - - VecView() = default; - friend class Vec; - friend class Vec::type>; - friend class VecView::type>; -}; - /* * Specialized vector implementation with multithreaded fill and uninitialized * memory optimizations. diff --git a/src/utilities/include/vec_view.h b/src/utilities/include/vec_view.h new file mode 100644 index 000000000..d4b4b649c --- /dev/null +++ b/src/utilities/include/vec_view.h @@ -0,0 +1,112 @@ +// Copyright 2023 The Manifold Authors. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include + +namespace manifold { + +/** + * View for Vec, can perform offset operation. + * This will be invalidated when the original vector is dropped or changes + * length. Roughly equivalent to std::span from c++20 + */ +template +class VecView { + public: + using Iter = T *; + using IterC = const T *; + + VecView(T *ptr_, int size_) : ptr_(ptr_), size_(size_) {} + + VecView(const VecView &other) { + ptr_ = other.ptr_; + size_ = other.size_; + } + + VecView &operator=(const VecView &other) { + ptr_ = other.ptr_; + size_ = other.size_; + return *this; + } + + // allows conversion to a const VecView + operator VecView() const { return {ptr_, size_}; } + + inline const T &operator[](int i) const { + if (i < 0 || i >= size_) throw std::out_of_range("Vec out of range"); + return ptr_[i]; + } + + inline T &operator[](int i) { + if (i < 0 || i >= size_) throw std::out_of_range("Vec out of range"); + return ptr_[i]; + } + + IterC cbegin() const { return ptr_; } + IterC cend() const { return ptr_ + size_; } + + IterC begin() const { return cbegin(); } + IterC end() const { return cend(); } + + Iter begin() { return ptr_; } + Iter end() { return ptr_ + size_; } + + const T &front() const { + if (size_ == 0) + throw std::out_of_range("attempt to take the front of an empty vector"); + return ptr_[0]; + } + + const T &back() const { + if (size_ == 0) + throw std::out_of_range("attempt to take the back of an empty vector"); + return ptr_[size_ - 1]; + } + + T &front() { + if (size_ == 0) + throw std::out_of_range("attempt to take the front of an empty vector"); + return ptr_[0]; + } + + T &back() { + if (size_ == 0) + throw std::out_of_range("attempt to take the back of an empty vector"); + return ptr_[size_ - 1]; + } + + int size() const { return size_; } + + bool empty() const { return size_ == 0; } + +#ifdef MANIFOLD_DEBUG + void Dump() { + std::cout << "Vec = " << std::endl; + for (int i = 0; i < size(); ++i) { + std::cout << i << ", " << ptr_[i] << ", " << std::endl; + } + std::cout << std::endl; + } +#endif + + protected: + T *ptr_ = nullptr; + int size_ = 0; + + VecView() = default; +}; + +} // namespace manifold \ No newline at end of file diff --git a/test/manifold_test.cpp b/test/manifold_test.cpp index 036f5182d..14884af8a 100644 --- a/test/manifold_test.cpp +++ b/test/manifold_test.cpp @@ -303,6 +303,23 @@ TEST(Manifold, Warp2) { EXPECT_NEAR(propBefore.volume, 321, 1); } +TEST(Manifold, WarpBatch) { + Manifold shape1 = + Manifold::Cube({2, 3, 4}).Warp([](glm::vec3& v) { v.x += v.z * v.z; }); + auto prop1 = shape1.GetProperties(); + + Manifold shape2 = + Manifold::Cube({2, 3, 4}).WarpBatch([](VecView vecs) { + for (glm::vec3& v : vecs) { + v.x += v.z * v.z; + } + }); + auto prop2 = shape2.GetProperties(); + + EXPECT_EQ(prop1.volume, prop2.volume); + EXPECT_EQ(prop1.surfaceArea, prop2.surfaceArea); +} + TEST(Manifold, Smooth) { Manifold tet = Manifold::Tetrahedron(); Manifold smooth = Manifold::Smooth(tet.GetMesh());