Add WarpBatch - array process warp with single callback

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()

bindings/python/manifold3d.cpp

@@ -70,38 +70,60 @@ namespace nb = nanobind;
 // helper to convert std::vector<glm::vec*> to numpy
 template <class T, int N, glm::qualifier Precision>
 nb::ndarray<nb::numpy, T, nb::shape<nb::any, N>> to_numpy(
- std::vector<glm::vec<N, T, Precision>> const &vecvec) {
+ glm::vec<N, T, Precision> const *begin,
+ glm::vec<N, T, Precision> 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<glm::vec*> to numpy
+template <class T, int N, glm::qualifier Precision>
+nb::ndarray<nb::numpy, T, nb::shape<nb::any, N>> to_numpy(
+ std::vector<glm::vec<N, T, Precision>> const &vec) {
+ return to_numpy(vec.data(), vec.data() + vec.size());
 }
 
 // helper to convert numpy to std::vector<glm::vec*>
 template <class T, size_t N, glm::qualifier Precision = glm::defaultp>
-std::vector<glm::vec<N, T, Precision>> to_glm_vector(
- nb::ndarray<nb::numpy, T, nb::shape<nb::any, N>> const &arr) {
- std::vector<glm::vec<N, T, Precision>> out;
- out.reserve(arr.shape(0));
+void to_glm_range(nb::ndarray<nb::numpy, T, nb::shape<nb::any, N>> const &arr,
+ glm::vec<int(N), T, Precision> *begin,
+ glm::vec<int(N), T, Precision> *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<glm::vec*>
+template <class T, size_t N, glm::qualifier Precision = glm::defaultp>
+void to_glm_vector(nb::ndarray<nb::numpy, T, nb::shape<nb::any, N>> const &arr,
+ std::vector<glm::vec<int(N), T, Precision>> &out) {
+ out.resize(arr.shape(0));
+ to_glm_range(arr, out.data(), out.data() + out.size());
 }
 
 using namespace manifold;
 
 typedef std::tuple<float, float> Float2;
 typedef std::tuple<float, float, float> Float3;
 
+using NumpyFloatNx2 = nb::ndarray<nb::numpy, float, nb::shape<nb::any, 2>>;
+using NumpyFloatNx3 = nb::ndarray<nb::numpy, float, nb::shape<nb::any, 3>>;
+using NumpyUintNx3 = nb::ndarray<nb::numpy, uint32_t, nb::shape<nb::any, 3>>;
+
 template <typename T>
 std::vector<T> toVector(const T *arr, size_t size) {
  return std::vector<T>(arr, arr + size);
@@ -150,12 +172,10 @@ NB_MODULE(manifold3d, m) {
  "triangulate",
  [](std::vector<nb::ndarray<nb::numpy, float, nb::shape<nb::any, 2>>>
  polys) {
- std::vector<std::vector<glm::vec2>> polys_vec;
- polys_vec.reserve(polys.size());
- for (auto &numpy : polys) {
- polys_vec.push_back(to_glm_vector(numpy));
+ std::vector<std::vector<glm::vec2>> 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<NumpyFloatNx3(NumpyFloatNx3)> &f) {
+ return self.WarpBatch([&f](VecView<glm::vec3> 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<NumpyFloatNx2(NumpyFloatNx2)> &f) {
+ return self.WarpBatch([&f](VecView<glm::vec2> 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 "

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<void(glm::vec2&)> warpFunc) const;
+ CrossSection WarpBatch(
+ std::function<void(VecView<glm::vec2>)> warpFunc) const;
  CrossSection Simplify(double epsilon = 1e-6) const;
 
  // Adapted from Clipper2 docs:

src/cross_section/src/cross_section.cpp

@@ -567,21 +567,43 @@ CrossSection CrossSection::Transform(const glm::mat3x2& m) const {
  */
 CrossSection CrossSection::Warp(
  std::function<void(glm::vec2&)> 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<glm::vec2> vecs) {
+ for (glm::vec2& p : vecs) {
+ warpFunc(p);
+ }
+ });
+}
+
+/**
+ * Same as CrossSection::Warp but calls warpFunc with begin and end
+ * iterators to all the vertices to be warped.
+ * Note that warpFunc begin/end pointers follow c++ iterator
+ * conventions - the end is exclusive and not to be dereferenced
+ *
+ * @param warpFunc A function that modifies multiple vertex positions.
+ */
+CrossSection CrossSection::WarpBatch(
+ std::function<void(VecView<glm::vec2>)> warpFunc) const {
+ std::vector<glm::vec2> tmp_verts;
+ C2::PathsD paths = GetPaths()->paths_; // deep copy
+ for (C2::PathD& path : paths) {
+ for (C2::PointD& p : path) {
+ tmp_verts.push_back(v2_of_pd(p));
+ }
+ }
+
+ warpFunc(VecView<glm::vec2>(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_)));
 }
 
 /**

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<void(glm::vec3&)>) const;
+ Manifold WarpBatch(std::function<void(VecView<glm::vec3>)>) const;
  Manifold SetProperties(
  int, std::function<void(float*, glm::vec3, const float*)>) const;
  Manifold CalculateCurvature(int gaussianIdx, int meanIdx) const;

src/manifold/src/impl.cpp

@@ -719,7 +719,14 @@ void Manifold::Impl::MarkFailure(Error status) {
 }
 
 void Manifold::Impl::Warp(std::function<void(glm::vec3&)> warpFunc) {
- thrust::for_each_n(thrust::host, vertPos_.begin(), NumVert(), warpFunc);
+ WarpBatch([&warpFunc](VecView<glm::vec3> vecs) {
+ thrust::for_each(thrust::host, vecs.begin(), vecs.end(), warpFunc);
+ });
+}
+
+void Manifold::Impl::WarpBatch(
+ std::function<void(VecView<glm::vec3>)> warpFunc) {
+ warpFunc(vertPos_.view());
  CalculateBBox();
  if (!IsFinite()) {
  MarkFailure(Error::NonFiniteVertex);

src/manifold/src/impl.h

@@ -76,6 +76,7 @@ struct Manifold::Impl {
  void Update();
  void MarkFailure(Error status);
  void Warp(std::function<void(glm::vec3&)> warpFunc);
+ void WarpBatch(std::function<void(VecView<glm::vec3>)> warpFunc);
  Impl Transform(const glm::mat4x3& transform) const;
  SparseIndices EdgeCollisions(const Impl& B, bool inverted = false) const;
  SparseIndices VertexCollisionsZ(VecView<const glm::vec3> vertsIn,

src/manifold/src/manifold.cpp

@@ -569,6 +569,20 @@ Manifold Manifold::Warp(std::function<void(glm::vec3&)> warpFunc) const {
  return Manifold(std::make_shared<CsgLeafNode>(pImpl));
 }
 
+/**
+ * Same as Manifold::Warp but calls warpFunc with begin and end
+ * iterators to all the vertices to be warped.
+ * Note that warpFunc begin/end pointers follow c++ iterator
+ * conventions - the end is exclusive and not to be dereferenced
+ * @param warpFunc A function that modifies multiple vertex positions.
+ */
+Manifold Manifold::WarpBatch(
+ std::function<void(VecView<glm::vec3>)> warpFunc) const {
+ auto pImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
+ pImpl->WarpBatch(warpFunc);
+ return Manifold(std::make_shared<CsgLeafNode>(pImpl));
+}
+
 /**
  * Create a new copy of this manifold with updated vertex properties by
  * supplying a function that takes the existing position and properties as