diff --git a/SRC/IsoOctree.h b/SRC/IsoOctree.h index 9a7cc37..a2089e0 100644 --- a/SRC/IsoOctree.h +++ b/SRC/IsoOctree.h @@ -158,6 +158,7 @@ class IsoOctree int setConforming(const std::vector& vertices,const std::vector >& polygons,const int& maxDepth,const int& setCenter,const Real& flatness,Point3D& translate,Real& scale,const int& noTransform); bool set(typename isoOctree::Octree::Traverser &traverser); + bool set(typename isoOctree::Octree::VectorizedTraverser &traverser); // A clean-up method to remove un-needed entries in the cornerValues hash-table void resetValues(void); diff --git a/SRC/api.cpp b/SRC/api.cpp index b59bc6a..45f9d3d 100644 --- a/SRC/api.cpp +++ b/SRC/api.cpp @@ -144,12 +144,9 @@ void PolygonToManifoldTriangleMesh( std::vector& vertices , const std::v } } } -} - -namespace isoOctree { -template -void buildMesh(typename isoOctree::Octree::Traverser &traverser, isoOctree::MeshInfo &output) +template +void buildMeshAny(Traverser &traverser, isoOctree::MeshInfo &output) { static bool isoOctreeCaseInit = false; if (!isoOctreeCaseInit) { @@ -189,11 +186,11 @@ void buildMesh(typename isoOctree::Octree::Traverser &traverser, isoOctree log_debug("Triangles: %zu", output.triangles.size()); } -template -void buildMeshWithPointCloudHint( - const std::function &)> &isoFunction, - const typename Octree::PointCloudHint &hint, - MeshInfo &output) +template +void buildMeshWithPointCloudHintAny( + const BuildMeshFunc &buildMeshFunc, + const typename isoOctree::Octree::PointCloudHint &hint, + isoOctree::MeshInfo &output) { log_debug("building mesh from point cloud hint with %zu point(s)", hint.nPoints); output.triangles.clear(); @@ -204,12 +201,15 @@ void buildMeshWithPointCloudHint( return; } - typename Octree::RootInfo root; + using Octree = isoOctree::Octree; + using Point3D = isoOctree::Point3D; + + typename Octree::RootInfo root; root.maxDepth = hint.maxDepth; { - Point3D minCorner = hint.points[0]; - Point3D maxCorner = hint.points[0]; + Point3D minCorner = hint.points[0]; + Point3D maxCorner = hint.points[0]; for (std::size_t i = 0; i < hint.nPoints; i++) { for (int j = 0; j < 3; j++) { @@ -236,7 +236,7 @@ void buildMeshWithPointCloudHint( }; } - using SearchTree = ZOrderOctree, Real>; + using SearchTree = ZOrderOctree; typename SearchTree::Parameters searchTreeParameters; // OK to cap. Just gets slower to search @@ -247,22 +247,14 @@ void buildMeshWithPointCloudHint( SearchTree searchTree(searchTreeParameters); searchTree.addData(hint.points, hint.nPoints); - class TraverserImpl : public Octree::Traverser { - private: - const std::function &)> &f; - const typename Octree::RootInfo &rootInfo; - const int subdivisionThreshold; - const SearchTree &searchTree; - - public: - const typename Octree::RootInfo &root() final { return rootInfo; } - - bool shouldExpand( - const typename Octree::Voxel &voxel, - const typename Octree::CornerValues &corners) final + buildMeshFunc( + root, + [&searchTree, &hint]( + const typename Octree::Voxel &voxel, + const typename Octree::CornerValues &corners) -> bool { (void)corners; - Point3D voxelCenter = { + Point3D voxelCenter = { Real(0.5) * voxel.width + voxel.minCorner[0], Real(0.5) * voxel.width + voxel.minCorner[1], Real(0.5) * voxel.width + voxel.minCorner[2] @@ -275,37 +267,143 @@ void buildMeshWithPointCloudHint( nMatches++; } - return nMatches >= subdivisionThreshold; + return nMatches >= hint.subdivisionThreshold; + }); +} +} + +namespace isoOctree { +template +void buildMesh(typename isoOctree::Octree::Traverser &traverser, isoOctree::MeshInfo &output) { + buildMeshAny(traverser, output); +} + +template +void buildMesh(typename isoOctree::Octree::VectorizedTraverser &traverser, isoOctree::MeshInfo &output) { + buildMeshAny(traverser, output); +} + +template +void buildMeshWithPointCloudHint( + const std::function &)> &isoFunction, + const typename Octree::PointCloudHint &hint, + MeshInfo &output) +{ + using Octree = Octree; + using RootInfo = typename Octree::RootInfo; + using IsoFunction = std::function &)>; + using ShouldExpandCallback = std::function; + + class TraverserImpl : public Octree::Traverser { + private: + const IsoFunction &isoFunction; + const ShouldExpandCallback &shouldExpandCallback; + const RootInfo &rootInfo; + + public: + const RootInfo &root() final { return rootInfo; } + + bool shouldExpand( + const typename Octree::Voxel &voxel, + const typename Octree::CornerValues &corners) final + { + return shouldExpandCallback(voxel, corners); } float isoValue(const Point3D &point) final { - return f(point); + return isoFunction(point); + } + + TraverserImpl( + const IsoFunction &isoFunction, + const ShouldExpandCallback &cb, + const RootInfo &r + ) : + isoFunction(isoFunction), + shouldExpandCallback(cb), + rootInfo(r) + {} + }; + + buildMeshWithPointCloudHintAny([&isoFunction, &output]( + const RootInfo &rootInfo, + const ShouldExpandCallback &shouldExpandCallback + ) { + TraverserImpl traverser(isoFunction, shouldExpandCallback, rootInfo); + buildMesh(traverser, output); + + }, hint, output); +} + +template +void buildMeshWithPointCloudHint( + const std::function> &, std::vector &)> &isoFunction, + const typename Octree::PointCloudHint &hint, + MeshInfo &output) +{ + using Octree = Octree; + using RootInfo = typename Octree::RootInfo; + using IsoFunction = std::function> &, std::vector &)>; + using ShouldExpandCallback = std::function; + + class TraverserImpl : public Octree::VectorizedTraverser { + private: + const IsoFunction &isoFunction; + const ShouldExpandCallback &shouldExpandCallback; + const RootInfo &rootInfo; + + public: + const RootInfo &root() final { return rootInfo; } + + bool shouldExpand( + const typename Octree::Voxel &voxel, + const typename Octree::CornerValues &corners) final + { + return shouldExpandCallback(voxel, corners); + } + + void isoValues(const std::vector> &points, std::vector &values) final { + isoFunction(points, values); } TraverserImpl( - const std::function &)> &f, - const typename Octree::RootInfo &r, - int subdivisionThreshold, - const SearchTree &searchTree + const IsoFunction &f, + const ShouldExpandCallback &cb, + const RootInfo &r ) : - f(f), - rootInfo(r), - subdivisionThreshold(subdivisionThreshold), - searchTree(searchTree) + isoFunction(f), + shouldExpandCallback(cb), + rootInfo(r) {} }; - TraverserImpl traverser(isoFunction, root, hint.subdivisionThreshold, searchTree); - buildMesh(traverser, output); + buildMeshWithPointCloudHintAny([&isoFunction, &output]( + const RootInfo &rootInfo, + const ShouldExpandCallback &shouldExpandCallback + ) { + TraverserImpl traverser(isoFunction, shouldExpandCallback, rootInfo); + buildMesh(traverser, output); + }, hint, output); } #define SPECIALIZE(real) \ template void buildMesh( \ Octree::Traverser &, \ MeshInfo &); \ +template void buildMesh( \ + Octree::VectorizedTraverser &, \ + MeshInfo &); \ template void buildMeshWithPointCloudHint( \ const std::function &)> &, \ const typename Octree::PointCloudHint &, \ + MeshInfo &); \ +template void buildMeshWithPointCloudHint( \ + const std::function> &, std::vector &)> &, \ + const typename Octree::PointCloudHint &, \ MeshInfo &) SPECIALIZE(float); @@ -327,6 +425,96 @@ bool IsoOctree::set(typename isoOctree::Octree return true; } +template +bool IsoOctree::set(typename isoOctree::Octree::VectorizedTraverser &traverser) { + maxDepth = traverser.root().maxDepth; + assert(maxDepth > 0); + cornerValues.clear(); + + std::vector keyList; + std::vector> evalList; + std::vector valueList; + + const auto &root = traverser.root(); + + struct Node { + typename OctNode::NodeIndex index; + OctNode* data; + }; + + std::vector nodes, nextNodes; + nodes.push_back({ + {}, // initialized to zero + &tree + }); + + for (int depth = 0; depth <= maxDepth; ++depth) { + keyList.clear(); + evalList.clear(); + + for (auto &node : nodes) { + const auto &nIdx = node.index; + + Point3D ctr; + Real w; + OctNode::CenterAndWidth(nIdx,ctr,w); + + for (int i=0; i::CornerIndex(nIdx, i, root.maxDepth); + if(cornerValues.find(key) == cornerValues.end()) { + int x,y,z; + Cube::FactorCornerIndex(i, x, y, z); + isoOctree::Point3D coords = getMappedCornerPosition( + root, + ctr.coords[0] + Real(x - 0.5) * w, + ctr.coords[1] + Real(y - 0.5) * w, + ctr.coords[2] + Real(z - 0.5) * w + ); + keyList.push_back(key); + evalList.push_back(coords); + cornerValues[key] = VertexData(); // will be evaluated on the next pass + } + } + } + + log_debug("level %d: evaluating %zu values", depth, keyList.size()); + valueList.resize(evalList.size()); + traverser.isoValues(evalList, valueList); + for (std::size_t i = 0; i < keyList.size(); ++i) { + cornerValues[keyList[i]] = VertexData(valueList[i]); + } + + if (depth == maxDepth) break; + + nextNodes.clear(); + for (auto &node : nodes) { + const auto &nIdx = node.index; + + typename isoOctree::Octree::CornerValues apiCornerVals; + for (int i=0; i::CornerIndex(nIdx, i, root.maxDepth); + apiCornerVals[i] = cornerValues.at(key).v; + } + + if (traverser.shouldExpand(convertIdx(nIdx, root), apiCornerVals)) { + if (!node.data->children) node.data->initChildren(); + for (int i=0; i::CornerIndex(nIdx, i, root.maxDepth); + + nextNodes.push_back({ + nIdx.child(i), + &node.data->children[i] + }); + } + } + } + + std::swap(nodes, nextNodes); + } + + return true; +} + template void IsoOctree::setChildren( OctNode* node, @@ -369,7 +557,6 @@ void IsoOctree::setChildren( if (!node->children) node->initChildren(); for (int i=0; i::CornerIndex(nIdx, i, root.maxDepth); setChildren(&node->children[i], nIdx.child(i), traverser); } } diff --git a/VERSION.txt b/VERSION.txt index 867e524..359a5b9 100644 --- a/VERSION.txt +++ b/VERSION.txt @@ -1 +1 @@ -1.2.0 \ No newline at end of file +2.0.0 \ No newline at end of file diff --git a/include/IsoOctree/IsoOctree.hpp b/include/IsoOctree/IsoOctree.hpp index 9190c09..b239565 100644 --- a/include/IsoOctree/IsoOctree.hpp +++ b/include/IsoOctree/IsoOctree.hpp @@ -75,16 +75,33 @@ template struct Octree { virtual float isoValue(const Point3D &point) = 0; // no virtual dtor, do not delete this class }; + + // evaluates the function one level at a time + struct VectorizedTraverser { + virtual const RootInfo &root() = 0; + virtual bool shouldExpand(const NodeIndex &node, const CornerValues &corners) = 0; + virtual void isoValues(const std::vector> &points, std::vector &result) = 0; + // no virtual dtor, do not delete this class + }; }; template ISO_OCTREE_API void buildMesh(typename Octree::Traverser &traverser, MeshInfo &output); +template +ISO_OCTREE_API void buildMesh(typename Octree::VectorizedTraverser &traverser, MeshInfo &output); + template ISO_OCTREE_API void buildMeshWithPointCloudHint( const std::function &)> &isoFunction, const typename Octree::PointCloudHint &hint, MeshInfo &output); + +template +ISO_OCTREE_API void buildMeshWithPointCloudHint( + const std::function> &, std::vector &)> &isoFunction, + const typename Octree::PointCloudHint &hint, + MeshInfo &output); } #endif diff --git a/python/bindings.cpp b/python/bindings.cpp index 9e6523d..e03c4e5 100644 --- a/python/bindings.cpp +++ b/python/bindings.cpp @@ -45,9 +45,33 @@ using RootInfo = Octree::RootInfo; using MeshInfo = isoOctree::MeshInfo; using ShouldExpandCallback = std::function; -using GetValueCallback = std::function; +using GetValueCallback = std::function; -class CallbackTraverser : public Octree::Traverser { +void convertCallback(GetValueCallback &pythonCallback, const std::vector &points, std::vector &values) { + const long n = static_cast(points.size()); + + py::array ret = pythonCallback(py::array_t( + std::vector { n, 3 }, + reinterpret_cast(points.data()) + )); + + if (ret.ndim() != 1) throw std::runtime_error("must return an 1-d array"); + if (ret.shape(0) != n) throw std::runtime_error("must the same number of values as there were input points"); + + ret = ret.attr("astype")(py::dtype("float32")); + + // Check if the array is already in row-major (C-contiguous) layout + if (!ret.attr("flags").attr("c_contiguous").cast()) { + // If not, make it row-major by copying the array + ret = ret.attr("copy")("C"); + } + + values.resize(n); + const float *data = reinterpret_cast(ret.data()); + std::copy(data, data + n, values.data()); +} + +class CallbackTraverser : public Octree::VectorizedTraverser { private: RootInfo rootInfo; ShouldExpandCallback shouldExpandCallback; @@ -69,8 +93,8 @@ class CallbackTraverser : public Octree::Traverser { return shouldExpandCallback(voxel, corners); } - float isoValue(const Point3D &point) final { - return getValueCallback(point); + void isoValues(const std::vector &points, std::vector &values) final { + convertCallback(getValueCallback, points, values); } }; } @@ -165,7 +189,15 @@ PYBIND11_MODULE(IsoOctree, m) { hint.subdivisionThreshold = subdivisionThreshold; std::unique_ptr> output = std::make_unique>(); - isoOctree::buildMeshWithPointCloudHint(isoFunction, hint, *output); + + std::function &, std::vector &)> cppCallback = [&isoFunction]( + const std::vector &points, + std::vector &values) + { + convertCallback(isoFunction, points, values); + }; + + isoOctree::buildMeshWithPointCloudHint(cppCallback, hint, *output); return std::move(output); }, py::arg("isoFunction"), diff --git a/python/examples/cayley.py b/python/examples/cayley.py index e4ddc8d..c709631 100644 --- a/python/examples/cayley.py +++ b/python/examples/cayley.py @@ -2,31 +2,31 @@ import numpy as np def getValue(p): - x, y, z = p + x, y, z = [p[:, i] for i in range(3)] # Cayley cubic return x**2 + y**2 - x**2*z + y**2*z + z**2 - 1 def getGradient(p): - x, y, z = p + x, y, z = [p[:, i] for i in range(3)] dx = 2*x - 2*x*z dy = 2*y + 2*y*z dz = -x**2 + y**2 + 2*z - return np.array([dx, dy, dz]) + return np.hstack([c[:, np.newaxis] for c in [dx, dy, dz]]) def getMinDot(node): h = node.width / 2 p = (node.minCorner[0] + h, node.minCorner[1] + h, node.minCorner[2] + h) normals = [] + grad_points = [] for dx in [-1, 1]: for dy in [-1, 1]: for dz in [-1, 1]: - p2 = (p[0] + dx * h, p[1] + dy * h, p[2] + dz * h) - grad = getGradient(p2) - normal = grad / max(np.linalg.norm(grad), 1e-6) - normals.append(normal) + grad_points.append([p[0] + dx * h, p[1] + dy * h, p[2] + dz * h]) + grad = getGradient(np.array(grad_points)) + normals = grad / np.maximum(np.linalg.norm(grad, axis=1), 1e-6)[:, np.newaxis] mean_normal = np.mean(normals, axis=0) mean_normal = mean_normal / max(np.linalg.norm(mean_normal), 1e-6) @@ -52,9 +52,8 @@ def shouldExpand(node, cornerValues): def refineVertices(vertices, iterations=5): for _ in range(iterations): - x, y, z = [vertices[:, i] for i in range(3)] - grad = getGradient((x, y, z)).T - value = getValue((x, y, z)) + grad = getGradient(vertices) + value = getValue(vertices) vertices = vertices - grad * ((value / np.sum(grad**2, axis=1)))[:, np.newaxis] return vertices diff --git a/python/examples/point_cloud.py b/python/examples/point_cloud.py index 020b25d..182a8e4 100644 --- a/python/examples/point_cloud.py +++ b/python/examples/point_cloud.py @@ -11,10 +11,8 @@ def buildMesh(points, normals, maxDepth, maxDist, kNearest, subdivisionThreshold tree = KDTree(points) md2 = maxDist**2 - def isoValue(point): - p0 = np.array(list(point)) + def isoValue(p0): _, ii = tree.query(p0, k=kNearest, distance_upper_bound=maxDist) - ii = [i for i in ii if i < len(points)] if len(ii) == 0: @@ -25,8 +23,11 @@ def isoValue(point): weights = 1.0 - np.array([min(d2, md2) / md2 for d2 in sqDists]) return np.sum([weights[i] * np.dot(normals[ii[i], :], p0 - points[ii[i], :]) for i in range(len(ii))]) + def isoValues(points): + return [isoValue(points[i, :]) for i in range(points.shape[0])] + return IsoOctree.buildMeshWithPointCloudHint( - isoValue, + isoValues, points, maxDepth=maxDepth, subdivisionThreshold=subdivisionThreshold) diff --git a/python/examples/simple.py b/python/examples/simple.py index 0f44eba..9ac8e8b 100644 --- a/python/examples/simple.py +++ b/python/examples/simple.py @@ -1,7 +1,7 @@ import IsoOctree def getValue(p): - x, y, z = p + x, y, z = [p[:, i] for i in range(3)] # Tanglecube return x*x*x*x - 5*x*x + y*y*y*y - 5*y*y + z*z*z*z - 5*z*z + 11.8