diff --git a/CMakeLists.txt b/CMakeLists.txt index 19f48136..d40966b7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -60,7 +60,7 @@ option(LIBIGL_CYCODEBASE "Build igl::cycodebase bindings" ON) FetchContent_Declare( libigl GIT_REPOSITORY https://github.com/libigl/libigl.git - GIT_TAG 72217e00bae6a7e041221e187aafbdc793f905d2 + GIT_TAG 0cbcaf886314ff331402c38fbb04cc4de8d7c033 ) FetchContent_MakeAvailable(libigl) diff --git a/include/default_types.h b/include/default_types.h index 85de9318..ddec29f9 100644 --- a/include/default_types.h +++ b/include/default_types.h @@ -1,3 +1,4 @@ +#pragma once #include #include #include diff --git a/include/parse_transforms.h b/include/parse_transforms.h new file mode 100644 index 00000000..d8307d1e --- /dev/null +++ b/include/parse_transforms.h @@ -0,0 +1,45 @@ +#pragma once +#include "default_types.h" +#include +#include +#include +#include + +namespace pyigl +{ + typedef Eigen::Transform AffineN3; + typedef std::vector > AffineN3List; + + /// Convert a Python list of 4×4 (or 3×4) matrices into the list of + /// Eigen::Transforms that the igl::swept_volume* functions expect. + /// + /// @param[in] transforms #transforms list of 4×4 or 3×4 matrices + /// @param[in] caller name used to prefix error messages + /// @return #transforms list of affine transformations + inline AffineN3List parse_transforms( + const std::vector &transforms, + const char * const caller) + { + if(transforms.empty()) + { + throw std::runtime_error( + std::string(caller)+": transforms must be non-empty"); + } + AffineN3List T; + T.reserve(transforms.size()); + for(const auto & M : transforms) + { + if(M.cols() != 4 || (M.rows() != 3 && M.rows() != 4)) + { + throw std::runtime_error( + std::string(caller)+": each transform must be 3×4 or 4×4"); + } + AffineN3 Ti; + // A 3×4 input leaves the implicit bottom row as [0 0 0 1] + Ti.matrix().setIdentity(); + Ti.matrix().topRows(M.rows()) = M; + T.emplace_back(Ti); + } + return T; + } +} diff --git a/pyproject.toml b/pyproject.toml index 58990358..146efdba 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,7 +13,7 @@ build-backend = "scikit_build_core.build" [project] name = "libigl" -version = "2.6.3.dev4" +version = "2.6.3.dev5" description = "libigl: A simple C++ geometry processing library" readme = "README.md" requires-python = ">=3.8" diff --git a/src/swept_volume.cpp b/src/swept_volume.cpp new file mode 100644 index 00000000..5e2c86a1 --- /dev/null +++ b/src/swept_volume.cpp @@ -0,0 +1,60 @@ +#include "default_types.h" +#include "parse_transforms.h" +#include +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + auto swept_volume( + const nb::DRef &V, + const nb::DRef &F, + /* Each transform is tiny; let nanobind copy them */ + const std::vector &transforms, + const igl::SignedDistanceType sign_type, + const Integer grid_res, + const Numeric isolevel) + { + const AffineN3List T = parse_transforms(transforms,"swept_volume"); + Eigen::MatrixXN SV; + Eigen::MatrixXI SF; + igl::swept_volume(V,F,T,sign_type,grid_res,isolevel,SV,SF); + return std::make_tuple(SV,SF); + } +} + +// Bind the wrapper to the Python module +void bind_swept_volume(nb::module_ &m) +{ + m.def( + "swept_volume", + &pyigl::swept_volume, + "V"_a, + "F"_a, + "transforms"_a, + "sign_type"_a = igl::SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER, + "grid_res"_a = 50, + "isolevel"_a = 0, +R"(Compute the surface of the swept volume of a solid object with surface +(V,F) mesh under going rigid motion. + +@param[in] V #V by 3 list of mesh positions in reference pose +@param[in] F #F by 3 list of mesh indices into rows of V +@param[in] transforms #transforms list of rigid transformations, one per time + step, each given as a 4×4 (or 3×4) matrix +@param[in] sign_type method for computing distance _sign_ +@param[in] grid_res number of grid cells on the longest side containing the + motion (enough cells to cover isolevel, plus one, will also be added on + each side as padding) +@param[in] isolevel distance level to be contoured as swept volume (in the + same units as V) +@return Tuple containing: + - SV: #SV by 3 list of mesh positions of the swept surface + - SF: #SF by 3 list of mesh faces into rows of SV)"); +} diff --git a/src/swept_volume_bounding_box.cpp b/src/swept_volume_bounding_box.cpp new file mode 100644 index 00000000..4fbfc7ba --- /dev/null +++ b/src/swept_volume_bounding_box.cpp @@ -0,0 +1,48 @@ +#include "default_types.h" +#include "parse_transforms.h" +#include +#include +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + auto swept_volume_bounding_box( + const nb::DRef &V, + /* Each transform is tiny; let nanobind copy them */ + const std::vector &transforms) + { + const AffineN3List T = + parse_transforms(transforms,"swept_volume_bounding_box"); + Eigen::AlignedBox box; + igl::swept_volume_bounding_box(V,T,box); + const Eigen::VectorXN min_corner = box.min(); + const Eigen::VectorXN max_corner = box.max(); + return std::make_tuple(min_corner,max_corner); + } +} + +// Bind the wrapper to the Python module +void bind_swept_volume_bounding_box(nb::module_ &m) +{ + m.def( + "swept_volume_bounding_box", + &pyigl::swept_volume_bounding_box, + "V"_a, + "transforms"_a, +R"(Construct an axis-aligned bounding box containing a shape undergoing a +motion sampled at a list of discrete rigid transformations. + +@param[in] V #V by 3 list of mesh positions in reference pose +@param[in] transforms #transforms list of rigid transformations, one per time + step, each given as a 4×4 (or 3×4) matrix +@return Tuple containing: + - min_corner: 3-vector of the minimum corner of the box + - max_corner: 3-vector of the maximum corner of the box)"); +} diff --git a/src/swept_volume_signed_distance.cpp b/src/swept_volume_signed_distance.cpp new file mode 100644 index 00000000..abb00939 --- /dev/null +++ b/src/swept_volume_signed_distance.cpp @@ -0,0 +1,91 @@ +#include "default_types.h" +#include "parse_transforms.h" +#include +#include +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + auto swept_volume_signed_distance( + const nb::DRef &V, + const nb::DRef &F, + /* Each transform is tiny; let nanobind copy them */ + const std::vector &transforms, + const igl::SignedDistanceType sign_type, + const nb::DRef &GV, + const nb::DRef &res, + const Numeric h, + const Numeric isolevel, + const nb::DRef &S0) + { + const AffineN3List T = + parse_transforms(transforms,"swept_volume_signed_distance"); + if(res.size() != 3) + { + throw std::runtime_error( + "swept_volume_signed_distance: res must be a 3-vector"); + } + if(res(0)*res(1)*res(2) != GV.rows()) + { + throw std::runtime_error( + "swept_volume_signed_distance: res(0)*res(1)*res(2) must equal GV.rows()"); + } + const Eigen::VectorXI r = res; + Eigen::VectorXN S; + if(S0.size() == 0) + { + igl::swept_volume_signed_distance(V,F,T,sign_type,GV,r,h,isolevel,S); + }else + { + if(S0.size() != GV.rows()) + { + throw std::runtime_error( + "swept_volume_signed_distance: S0 must have GV.rows() entries"); + } + igl::swept_volume_signed_distance(V,F,T,sign_type,GV,r,h,isolevel,S0,S); + } + return S; + } +} + +// Bind the wrapper to the Python module +void bind_swept_volume_signed_distance(nb::module_ &m) +{ + m.def( + "swept_volume_signed_distance", + &pyigl::swept_volume_signed_distance, + "V"_a, + "F"_a, + "transforms"_a, + "sign_type"_a = igl::SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER, + "GV"_a, + "res"_a, + "h"_a, + "isolevel"_a = std::numeric_limits::infinity(), + "S0"_a = Eigen::VectorXN(), +R"(Compute the signed distance to a sweep surface of a mesh under-going a +rigid motion discretely sampled at a list of transformations at a grid. + +@param[in] V #V by 3 list of mesh positions in reference pose +@param[in] F #F by 3 list of triangle indices [0,n) +@param[in] transforms #transforms list of rigid transformations, one per time + step, each given as a 4×4 (or 3×4) matrix +@param[in] sign_type method for computing distance _sign_ +@param[in] GV #GV by 3 list of evaluation point grid positions +@param[in] res 3-long resolution of the GV grid +@param[in] h edge-length of grid +@param[in] isolevel isolevel to "focus" on; grid positions far enough away + from isolevel (based on h) will get approximate values. Set + isolevel=inf (the default) to get good values everywhere (slow and + unnecessary if just trying to extract the isolevel-level set). +@param[in] S0 #GV list of initial values (the minimum with these is taken); + empty (the default) to start from scratch +@return S #GV list of signed distances)"); +} diff --git a/src/voxel_grid.cpp b/src/voxel_grid.cpp index 934166aa..0bb247f3 100644 --- a/src/voxel_grid.cpp +++ b/src/voxel_grid.cpp @@ -3,6 +3,8 @@ #include #include #include +#include +#include namespace nb = nanobind; using namespace nb::literals; @@ -20,6 +22,28 @@ namespace pyigl igl::voxel_grid(V, offset, s, pad_count, GV, side); return std::make_tuple(GV, side); } + + auto voxel_grid_box( + const nb::DRef &min_corner, + const nb::DRef &max_corner, + const int s, + const int pad_count) + { + if(min_corner.size() != 3 || max_corner.size() != 3) + { + throw std::runtime_error( + "voxel_grid: min_corner and max_corner must be 3-vectors"); + } + // Named locals: `AlignedBox box(Matrix(a),Matrix(b));` would parse as a + // function declaration + const Eigen::Matrix mn = min_corner; + const Eigen::Matrix mx = max_corner; + const Eigen::AlignedBox box(mn,mx); + Eigen::MatrixXN GV; + Eigen::VectorXI side; + igl::voxel_grid(box, s, pad_count, GV, side); + return std::make_tuple(GV, side); + } } // Bind the wrappers to the Python module @@ -39,4 +63,19 @@ void bind_voxel_grid(nb::module_ &m) @param[in] s Number of cell centers on the largest side @param[in] pad_count Number of cells beyond the box @return Tuple (GV, side) where GV contains cell center positions and side defines grid dimensions)"); + m.def( + "voxel_grid", + &pyigl::voxel_grid_box, + "min_corner"_a, + "max_corner"_a, + "s"_a, + "pad_count"_a=0, + R"(Construct the cell center positions of a regular voxel grid (lattice) +made of perfectly square voxels enclosing a given axis-aligned box. + +@param[in] min_corner 3-vector of the minimum corner of the box to enclose +@param[in] max_corner 3-vector of the maximum corner of the box to enclose +@param[in] s Number of cell centers on the largest side (including 2*pad_count) +@param[in] pad_count Number of cells beyond the box +@return Tuple (GV, side) where GV contains cell center positions and side defines grid dimensions)"); } diff --git a/tests/test_all.py b/tests/test_all.py index ef5e051c..a3f5ba7d 100644 --- a/tests/test_all.py +++ b/tests/test_all.py @@ -359,6 +359,25 @@ def test_voxel(): GV,side = igl.voxel_grid(V,s=10) GV,side = igl.voxel_grid(V,s=10,offset=0.1,pad_count=2) +def test_voxel_grid_box(): + V,_,_ = single_tet() + min_corner = V.min(axis=0) + max_corner = V.max(axis=0) + GV,side = igl.voxel_grid(min_corner,max_corner,s=10,pad_count=2) + assert GV.dtype == np.float64 + assert side.dtype == np.int64 + assert GV.shape == (np.prod(side),3) + assert side.shape == (3,) + # Enclosing the corners is the same as enclosing the points themselves + GV2,side2 = igl.voxel_grid(V,0.0,s=10,pad_count=2) + assert np.allclose(GV,GV2) + assert np.array_equal(side,side2) + # Positional args resolve to the box overload too + GV3,side3 = igl.voxel_grid(min_corner,max_corner,10,2) + assert np.allclose(GV,GV3) + with pytest.raises(RuntimeError): + igl.voxel_grid(min_corner[:2],max_corner,s=10) + def test_sample(): V,F = igl.icosahedron() @@ -1958,3 +1977,92 @@ def test_vertex_components_from_adjacency_matrix(): # counts is per-component and sums to the number of vertices assert counts.sum() == n assert sorted(counts.ravel().tolist()) == [2, 3] + +def test_swept_volume(): + V,F = igl.icosahedron() + # Rigid motion: translate 2 units along x in 5 steps + T = np.tile(np.eye(4,dtype=np.float64),(5,1,1)) + T[:,0,3] = np.linspace(0,2,5) + SV,SF = igl.swept_volume(V,F,T,igl.SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER,16,0.0) + assert SV.dtype == np.float64 + assert SF.dtype == np.int64 + assert SV.shape[1] == 3 + assert SF.shape[1] == 3 + assert SV.shape[0] > 0 + assert SF.shape[0] > 0 + # Swept volume spans the motion (icosahedron has unit-ish radius) + assert SV[:,0].min() < -0.5 + assert SV[:,0].max() > 2.5 + # A list of 3×4 transforms is accepted and equivalent + SV2,SF2 = igl.swept_volume(V,F,[Ti[:3,:] for Ti in T],grid_res=16) + assert np.allclose(SV,SV2) + assert np.array_equal(SF,SF2) + # A positive isolevel dilates the result + SV3,SF3 = igl.swept_volume(V,F,T,grid_res=16,isolevel=0.25) + assert SV3[:,1].max() > SV[:,1].max() + with pytest.raises(RuntimeError): + igl.swept_volume(V,F,[]) + with pytest.raises(RuntimeError): + igl.swept_volume(V,F,[np.eye(3)]) + +def swept_translation(n=5,dist=2.0): + """#n by 4 by 4 list of transforms translating `dist` along x.""" + T = np.tile(np.eye(4,dtype=np.float64),(n,1,1)) + T[:,0,3] = np.linspace(0,dist,n) + return T + +def test_swept_volume_bounding_box(): + V,F = igl.icosahedron() + T = swept_translation() + mn,mx = igl.swept_volume_bounding_box(V,T) + assert mn.dtype == np.float64 + assert mx.dtype == np.float64 + assert mn.shape == (3,) + assert mx.shape == (3,) + assert np.all(mn <= mx) + # Box is the reference pose's box swept 2 units along x + assert np.allclose(mn,V.min(axis=0)) + assert np.allclose(mx,V.max(axis=0)+np.array([2.0,0,0])) + # 3x4 transforms give the same box + mn2,mx2 = igl.swept_volume_bounding_box(V,[Ti[:3,:] for Ti in T]) + assert np.allclose(mn,mn2) + assert np.allclose(mx,mx2) + with pytest.raises(RuntimeError): + igl.swept_volume_bounding_box(V,[]) + +def test_swept_volume_signed_distance(): + V,F = igl.icosahedron() + T = swept_translation() + grid_res = 16 + isolevel = 0.0 + mn,mx = igl.swept_volume_bounding_box(V,T) + h = (mx-mn).max()/(grid_res-1) + pad = max(int(np.ceil(isolevel/h)),0)+1 + GV,res = igl.voxel_grid(mn,mx,s=grid_res+2*pad,pad_count=pad) + S = igl.swept_volume_signed_distance( + V,F,T,igl.SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER,GV,res,h,isolevel) + assert S.dtype == np.float64 + assert S.shape == (GV.shape[0],) + assert not np.any(np.isnan(S)) + # Signs straddle the surface + assert S.min() < 0 and S.max() > 0 + # Contouring these values reproduces igl.swept_volume exactly + SV,SF,_ = igl.marching_cubes( + S-isolevel,GV,int(res[0]),int(res[1]),int(res[2]),0.0) + SV2,SF2 = igl.swept_volume( + V,F,T,igl.SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER,grid_res,isolevel) + assert np.allclose(SV,SV2) + assert np.array_equal(SF,SF2) + # Feeding the result back in as S0 takes a min with itself: idempotent + S2 = igl.swept_volume_signed_distance( + V,F,T,igl.SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER,GV,res,h,isolevel,S) + assert np.allclose(S,S2) + # isolevel=inf (the default) is exact everywhere, so no NaNs to flood fill + Sinf = igl.swept_volume_signed_distance(V,F,T,GV=GV,res=res,h=h) + assert not np.any(np.isnan(Sinf)) + near = np.abs(Sinf) < h + assert np.allclose(S[near],Sinf[near]) + with pytest.raises(RuntimeError): + igl.swept_volume_signed_distance(V,F,T,GV=GV,res=res[:2],h=h) + with pytest.raises(RuntimeError): + igl.swept_volume_signed_distance(V,F,T,GV=GV,res=res,h=h,S0=S[:5])