diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..aeeac98 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,68 @@ +name: CI + +on: + push: + branches: ["**"] + pull_request: + branches: ["**"] +jobs: + tester: + name: Testing + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, macos-latest, macos-14,windows-latest] + python-version: ${{ fromJSON(vars.PYTHON_VERSIONS)}} + env: + MACOSX_DEPLOYMENT_TARGET: '15.0' + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + - if: matrix.os == 'ubuntu-latest' + run: | + sudo apt-get install -y libgl1-mesa-dev libeigen3-dev libcgal-dev + - if: ${{ startsWith(matrix.os, 'macos') }} + run: | + brew install cgal + # Cache vcpkg on Windows + - name: Restore vcpkg cache + if: runner.os == 'Windows' + id: cache-vcpkg + uses: actions/cache@v4 + with: + path: | + C:/vcpkg + key: vcpkg-${{ runner.os }}-${{ hashFiles('vcpkg.json') }} + restore-keys: | + vcpkg-${{ runner.os }}- + - name: Install vcpkg dependencies (Windows only) + if: runner.os == 'Windows' + shell: pwsh + run: | + $Env:VCPKG_ROOT = "C:/vcpkg" + if (!(Test-Path $Env:VCPKG_ROOT)) { + git clone https://github.com/microsoft/vcpkg $Env:VCPKG_ROOT + & "$Env:VCPKG_ROOT/bootstrap-vcpkg.bat" + } + # Only install if packages are missing + if (!(Test-Path "$Env:VCPKG_ROOT/installed/x64-windows/cgal")) { + & "$Env:VCPKG_ROOT/vcpkg.exe" install yasm-tool cgal + } + echo "VCPKG_ROOT=$Env:VCPKG_ROOT" >> $Env:GITHUB_ENV + # now inject your extra CMake flags: + $toolchain = "$Env:VCPKG_ROOT/scripts/buildsystems/vcpkg.cmake" + echo "CMAKE_ARGS=-DCMAKE_TOOLCHAIN_FILE=${toolchain} -DVCPKG_ROOT=${Env:VCPKG_ROOT}" >> $Env:GITHUB_ENV + echo $CMAKE_ARGS + - name: Pip install + run: | + pip install .[test] + - name: Run tests + run: | + pytest + diff --git a/.gitignore b/.gitignore index 2d5550d..225d4d8 100644 --- a/.gitignore +++ b/.gitignore @@ -43,4 +43,5 @@ Thumbs.db loop_cgal/__pycache__/__init__.cpython-311.pyc build/* *.vtk -*.ply \ No newline at end of file +*.ply +build_win64/* \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 824be76..fe9c2a8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,26 +25,30 @@ add_library(_loop_cgal MODULE target_link_libraries(_loop_cgal PRIVATE pybind11::module CGAL::CGAL) target_include_directories(_loop_cgal PRIVATE ${CMAKE_SOURCE_DIR}/src) -set_target_properties(_loop_cgal PROPERTIES PREFIX "" SUFFIX ".pyd") +if(WIN32) + set_target_properties(_loop_cgal PROPERTIES PREFIX "" SUFFIX ".pyd") +else() + set_target_properties(_loop_cgal PROPERTIES PREFIX "" SUFFIX ".so") +endif() # Windows: copy required DLLs from VCPKG_ROOT + + if(WIN32) - if(NOT DEFINED ENV{VCPKG_ROOT}) - message(FATAL_ERROR "VCPKG_ROOT environment variable is not set!") - endif() - set(VCPKG_ROOT $ENV{VCPKG_ROOT}) - message(STATUS "Using VCPKG_ROOT: ${VCPKG_ROOT}") - - add_custom_command(TARGET _loop_cgal POST_BUILD - COMMAND ${CMAKE_COMMAND} -E copy_if_different - "${VCPKG_ROOT}/packages/gmp_x64-windows/bin/gmp-10.dll" - "$" - COMMAND ${CMAKE_COMMAND} -E copy_if_different - "${VCPKG_ROOT}/packages/mpfr_x64-windows/bin/mpfr-6.dll" - "$" - ) + install(CODE " + file(INSTALL + DESTINATION \"\${CMAKE_INSTALL_PREFIX}/loop_cgal\" + TYPE FILE + FILES + \"${VCPKG_ROOT}/packages/gmp_x64-windows/bin/gmp-10.dll\" + \"${VCPKG_ROOT}/packages/mpfr_x64-windows/bin/mpfr-6.dll\" + ) + message(\"[CMake Install] Copied gmp-10.dll and mpfr-6.dll to install directory.\") + ") endif() + + # Install for Python packages install(TARGETS _loop_cgal LIBRARY DESTINATION ${CMAKE_INSTALL_PREFIX}/loop_cgal diff --git a/examples/clip_test.py b/examples/clip_example.py similarity index 100% rename from examples/clip_test.py rename to examples/clip_example.py diff --git a/examples/cut_example.py b/examples/cut_example.py new file mode 100644 index 0000000..cf6ebb3 --- /dev/null +++ b/examples/cut_example.py @@ -0,0 +1,34 @@ +import numpy as np +from loop_cgal import TriMesh, set_verbose +from LoopStructural.datatypes import BoundingBox +set_verbose(True) + + +bb = BoundingBox(np.zeros(3), np.ones(3)) +grid = bb.structured_grid().vtk() +grid["scalars"] = grid.points[:, 0] +surface = grid.contour([0.5]) +property = surface.points[:, 1].copy() +surface['prop'] = property +surface.save("input_mesh.vtk") + +mesh = TriMesh(surface) + +# Define a scalar property that is the Y coordinate of each vertex + +# Cut value at y = 0.6; this will split triangles crossing y=0.6 +cut_value = property.mean() +# mesh.remesh(0.10, protect_constraints=False) + +mesh.vtk().save('before_cut.vtk') +# Invoke the new method (property is python list or numpy array) +mesh.cut_with_implicit_function(property.tolist(), cut_value) +mesh.reverse_face_orientation() +# Convert back to pyvista for visualization +out = mesh.to_pyvista() +out.save('cut_mesh.vtk') +# # Visualize the result +# pl = pv.Plotter() +# pl.add_mesh(out, show_edges=True, color='lightgray') +# pl.add_mesh(pv.PolyData(np.array([[cut_value, -1e3, 0.0],[cut_value,1e3,0.0]])), color='red') +# pl.show() diff --git a/examples/split_test.py b/examples/split_example.py similarity index 100% rename from examples/split_test.py rename to examples/split_example.py diff --git a/loop_cgal/__init__.py b/loop_cgal/__init__.py index a52d5e6..fcbbe20 100644 --- a/loop_cgal/__init__.py +++ b/loop_cgal/__init__.py @@ -7,104 +7,86 @@ import pyvista as pv from ._loop_cgal import TriMesh as _TriMesh -from ._loop_cgal import verbose # noqa: F401 +from ._loop_cgal import verbose # noqa: F401 from ._loop_cgal import set_verbose as set_verbose - - - - -def validate_pyvista_polydata(surface: pv.PolyData, surface_name: str = "surface") -> None: - """Validate a PyVista PolyData object. - - Parameters - ---------- - surface : pv.PolyData - The surface to validate - surface_name : str - Name of the surface for error messages - - Raises - ------ - ValueError - If the surface is invalid - """ - if not isinstance(surface, pv.PolyData): - raise ValueError(f"{surface_name} must be a pyvista.PolyData object") - - if surface.n_points == 0: - raise ValueError(f"{surface_name} has no points") - - if surface.n_cells == 0: - raise ValueError(f"{surface_name} has no cells") - - points = np.asarray(surface.points) - if not np.isfinite(points).all(): - raise ValueError(f"{surface_name} points contain NaN or infinite values") - +from .utils import validate_pyvista_polydata, validate_vertices_and_faces class TriMesh(_TriMesh): """ A class for handling triangular meshes using CGAL. - + Inherits from the base TriMesh class and provides additional functionality. """ + def __init__(self, surface: pv.PolyData): # Validate input surface validate_pyvista_polydata(surface, "input surface") - + # Triangulate to ensure we have triangular faces surface = surface.triangulate() - + # Extract vertices and triangles verts = np.array(surface.points, dtype=np.float64).copy() faces = surface.faces.reshape(-1, 4)[:, 1:].copy().astype(np.int32) - - # Additional validation on extracted data - if verts.size == 0: - raise ValueError("Surface has no vertices after triangulation") - - if faces.size == 0: - raise ValueError("Surface has no triangular faces after triangulation") - - if not np.isfinite(verts).all(): - raise ValueError("Surface vertices contain NaN or infinite values") - - # Check triangle indices - max_vertex_index = verts.shape[0] - 1 - if faces.min() < 0: - raise ValueError("Surface has negative triangle indices") - - if faces.max() > max_vertex_index: - raise ValueError(f"Surface triangle indices exceed vertex count (max index: {faces.max()}, vertex count: {verts.shape[0]})") - # Check for degenerate triangles - # build a ntris x nverts matrix - # populate with true for vertex in each triangle - # sum rows and if not equal to 3 then it is degenerate - face_idx = np.arange(faces.shape[0]) - face_idx = np.tile(face_idx, (3,1)).T.flatten() - faces_flat = faces.flatten() - m = sp.coo_matrix( - (np.ones(faces_flat.shape[0]), (faces_flat, face_idx)), - shape=(verts.shape[0], faces.shape[0]), - dtype=bool, - ) - # coo duplicates entries so just make sure its boolean - m = m > 0 - if not np.all(m.sum(axis=0) == 3): - degen_idx = np.where(m.sum(axis=0) != 3)[1] - raise ValueError(f"Surface contains degenerate triangles: {degen_idx} (each triangle must have exactly 3 vertices)") - + if not validate_vertices_and_faces(verts, faces): + raise ValueError("Invalid surface geometry") super().__init__(verts, faces) - - def to_pyvista(self, area_threshold: float = 1e-6, # this is the area threshold for the faces, if the area is smaller than this it will be removed - duplicate_vertex_threshold: float = 1e-4, # this is the threshold for duplicate vertices - ) -> pv.PolyData: + @classmethod + def from_vertices_and_triangles( + cls, vertices: np.ndarray, triangles: np.ndarray + ) -> TriMesh: + """ + Create a TriMesh from vertices and triangle indices. + + Parameters + ---------- + vertices : np.ndarray + An array of shape (n_vertices, 3) containing the vertex coordinates. + triangles : np.ndarray + An array of shape (n_triangles, 3) containing the triangle vertex indices. + + Returns + ------- + TriMesh + The created TriMesh object. + """ + # Create a temporary PyVista PolyData object for validation + if not validate_vertices_and_faces(vertices, triangles): + raise ValueError("Invalid vertices or triangles") + surface = pv.PolyData(vertices, np.hstack((np.full((triangles.shape[0], 1), 3), triangles)).flatten()) + return cls(surface) + + def get_vertices_and_triangles( + self, + area_threshold: float = 1e-6, # this is the area threshold for the faces, if the area is smaller than this it will be removed + duplicate_vertex_threshold: float = 1e-4, + ) -> Tuple[np.ndarray, np.ndarray]: + """ + Get the vertices and triangle indices of the TriMesh. + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + A tuple containing: + - An array of shape (n_vertices, 3) with the vertex coordinates. + - An array of shape (n_triangles, 3) with the triangle vertex indices. + """ + np_mesh = self.save(area_threshold, duplicate_vertex_threshold) + vertices = np.array(np_mesh.vertices).copy() + triangles = np.array(np_mesh.triangles).copy() + return vertices, triangles + + def to_pyvista( + self, + area_threshold: float = 1e-6, # this is the area threshold for the faces, if the area is smaller than this it will be removed + duplicate_vertex_threshold: float = 1e-4, # this is the threshold for duplicate vertices + ) -> pv.PolyData: """ Convert the TriMesh to a pyvista PolyData object. - + Returns ------- pyvista.PolyData @@ -115,3 +97,23 @@ def to_pyvista(self, area_threshold: float = 1e-6, # this is the area threshold triangles = np.array(np_mesh.triangles).copy() return pv.PolyData.from_regular_faces(vertices, triangles) + def vtk( + self, + area_threshold: float = 1e-6, + duplicate_vertex_threshold: float = 1e-4, + ) -> pv.PolyData: + """ + Alias for to_pyvista method. + """ + return self.to_pyvista(area_threshold, duplicate_vertex_threshold) + + def copy(self) -> TriMesh: + """ + Create a copy of the TriMesh. + + Returns + ------- + TriMesh + A copy of the TriMesh object. + """ + return TriMesh(self.to_pyvista()) diff --git a/loop_cgal/bindings.cpp b/loop_cgal/bindings.cpp index 72fc1c8..1546f4b 100644 --- a/loop_cgal/bindings.cpp +++ b/loop_cgal/bindings.cpp @@ -14,12 +14,18 @@ PYBIND11_MODULE(_loop_cgal, m) .def(py::init<>()) .def_readwrite("vertices", &NumpyMesh::vertices) .def_readwrite("triangles", &NumpyMesh::triangles); + py::enum_(m, "ImplicitCutMode") + .value("PRESERVE_INTERSECTION", ImplicitCutMode::PRESERVE_INTERSECTION) + .value("KEEP_POSITIVE_SIDE", ImplicitCutMode::KEEP_POSITIVE_SIDE) + .value("KEEP_NEGATIVE_SIDE", ImplicitCutMode::KEEP_NEGATIVE_SIDE) + .export_values(); py::class_(m, "TriMesh") .def(py::init &, const pybind11::array_t &>(), py::arg("vertices"), py::arg("triangles")) .def("cut_with_surface", &TriMesh::cutWithSurface, py::arg("surface"), py::arg("preserve_intersection") = false, - py::arg("preserve_intersection_clipper") = false) + py::arg("preserve_intersection_clipper") = false, + py::arg("use_exact_kernel") = true) .def("remesh", &TriMesh::remesh, py::arg("split_long_edges") = true, py::arg("target_edge_length") = 10.0, py::arg("number_of_iterations") = 3, @@ -31,6 +37,9 @@ PYBIND11_MODULE(_loop_cgal, m) "Reverse the face orientation of the mesh.") .def("add_fixed_edges", &TriMesh::add_fixed_edges, py::arg("pairs"), - "Vertex index pairs defining edges to be fixed in mesh when remeshing."); + "Vertex index pairs defining edges to be fixed in mesh when remeshing.") + .def("cut_with_implicit_function", &TriMesh::cut_with_implicit_function, + py::arg("property"), py::arg("value"),py::arg("cutmode") = ImplicitCutMode::KEEP_POSITIVE_SIDE, + "Cut the mesh with an implicit function defined by vertex properties."); } // End of PYBIND11_MODULE \ No newline at end of file diff --git a/loop_cgal/utils.py b/loop_cgal/utils.py new file mode 100644 index 0000000..011395c --- /dev/null +++ b/loop_cgal/utils.py @@ -0,0 +1,98 @@ +import pyvista as pv +import numpy as np +import scipy.sparse as sp +def validate_pyvista_polydata( + surface: pv.PolyData, surface_name: str = "surface" +) -> None: + """Validate a PyVista PolyData object. + + Parameters + ---------- + surface : pv.PolyData + The surface to validate + surface_name : str + Name of the surface for error messages + + Raises + ------ + ValueError + If the surface is invalid + """ + if not isinstance(surface, pv.PolyData): + raise ValueError(f"{surface_name} must be a pyvista.PolyData object") + + if surface.n_points == 0: + raise ValueError(f"{surface_name} has no points") + + if surface.n_cells == 0: + raise ValueError(f"{surface_name} has no cells") + + points = np.asarray(surface.points) + if not np.isfinite(points).all(): + raise ValueError(f"{surface_name} points contain NaN or infinite values") + + +def validate_vertices_and_faces(verts, faces): + """Validate vertices and faces arrays. + Parameters + ---------- + verts : np.ndarray + + An array of shape (n_vertices, 3) containing the vertex coordinates. + faces : np.ndarray + + An array of shape (n_faces, 3) containing the triangle vertex indices. + Raises + ------ + ValueError + If the vertices or faces are invalid. + """ + if type(verts) is not np.ndarray: + try: + verts = np.array(verts) + except Exception: + raise ValueError("Vertices must be a numpy array") + if type(faces) is not np.ndarray: + try: + faces = np.array(faces) + except Exception: + raise ValueError("Faces must be a numpy array") + # Additional validation on extracted data + if verts.size == 0: + raise ValueError("Surface has no vertices after triangulation") + + if faces.size == 0: + raise ValueError("Surface has no triangular faces after triangulation") + + if not np.isfinite(verts).all(): + raise ValueError("Surface vertices contain NaN or infinite values") + + # Check triangle indices + max_vertex_index = verts.shape[0] - 1 + if faces.min() < 0: + raise ValueError("Surface has negative triangle indices") + + if faces.max() > max_vertex_index: + raise ValueError( + f"Surface triangle indices exceed vertex count (max index: {faces.max()}, vertex count: {verts.shape[0]})" + ) + # Check for degenerate triangles + # build a ntris x nverts matrix + # populate with true for vertex in each triangle + # sum rows and if not equal to 3 then it is degenerate + face_idx = np.arange(faces.shape[0]) + face_idx = np.tile(face_idx, (3, 1)).T.flatten() + faces_flat = faces.flatten() + m = sp.coo_matrix( + (np.ones(faces_flat.shape[0]), (faces_flat, face_idx)), + shape=(verts.shape[0], faces.shape[0]), + dtype=bool, + ) + # coo duplicates entries so just make sure its boolean + m = m > 0 + if not np.all(m.sum(axis=0) == 3): + degen_idx = np.where(m.sum(axis=0) != 3)[1] + raise ValueError( + f"Surface contains degenerate triangles: {degen_idx} (each triangle must have exactly 3 vertices)" + ) + return True \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 1cb355b..21daf0f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,7 +26,10 @@ classifiers = [ "Programming Language :: Python :: 3.12", ] -dependencies = ['pyvista','scipy','numpy'] +dependencies = ['pyvista','vtk','scipy','numpy'] + +[project.optional-dependencies] +test = ['pytest'] [tool.scikit-build] wheel.expand-macos-universal-tags = true @@ -83,3 +86,5 @@ manylinux-x86_64-image = "lachlangrose/manylinuxeigencgal" [tool.cibuildwheel.macos] archs = ["x86_64","arm64"] +[tool.setuptools.package-data] +"loop_cgal" = ["*.dll"] \ No newline at end of file diff --git a/src/mesh.cpp b/src/mesh.cpp index d679a25..157b2f3 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -63,14 +63,16 @@ TriMesh::TriMesh(const pybind11::array_t &vertices, _mesh.add_vertex(Point(verts(i, 0), verts(i, 1), verts(i, 2)))); } - for (ssize_t i = 0; i < tris.shape(0); ++i) { + for (ssize_t i = 0; i < tris.shape(0); ++i) + { int v0 = tris(i, 0); int v1 = tris(i, 1); int v2 = tris(i, 2); // Check that all vertex indices are valid if (v0 < 0 || v0 >= vertex_indices.size() || v1 < 0 || - v1 >= vertex_indices.size() || v2 < 0 || v2 >= vertex_indices.size()) { + v1 >= vertex_indices.size() || v2 < 0 || v2 >= vertex_indices.size()) + { std::cerr << "Warning: Triangle " << i << " has invalid vertex indices: (" << v0 << ", " << v1 << ", " << v2 << "). Skipping." << std::endl; @@ -78,7 +80,8 @@ TriMesh::TriMesh(const pybind11::array_t &vertices, } // Check for degenerate triangles - if (v0 == v1 || v1 == v2 || v0 == v2) { + if (v0 == v1 || v1 == v2 || v0 == v2) + { std::cerr << "Warning: Triangle " << i << " is degenerate: (" << v0 << ", " << v1 << ", " << v2 << "). Skipping." << std::endl; continue; @@ -105,7 +108,6 @@ void TriMesh::init() { _fixedEdges = collect_border_edges(_mesh); - if (LoopCGAL::verbose) { std::cout << "Found " << _fixedEdges.size() << " fixed edges." << std::endl; @@ -141,7 +143,7 @@ void TriMesh::add_fixed_edges(const pybind11::array_t &pairs) << std::endl; continue; } - if (!_mesh.is_valid(edge)) // Check if the halfedge is valid + if (!_mesh.is_valid(edge)) // Check if the halfedge is valid { std::cerr << "Invalid half-edge for vertices (" << v0 << ", " << v1 << ")" << std::endl; @@ -272,42 +274,54 @@ void TriMesh::remesh(bool split_long_edges, void TriMesh::reverseFaceOrientation() { // Reverse the face orientation of the mesh + PMP::remove_isolated_vertices(_mesh); + if (!CGAL::is_valid_polygon_mesh(_mesh, LoopCGAL::verbose)) + { + std::cerr << "Mesh is not valid before reversing face orientations." + << std::endl; + return; + } PMP::reverse_face_orientations(_mesh); if (!CGAL::is_valid_polygon_mesh(_mesh, LoopCGAL::verbose)) { std::cerr << "Mesh is not valid after reversing face orientations." << std::endl; } - } void TriMesh::cutWithSurface(TriMesh &clipper, bool preserve_intersection, - bool preserve_intersection_clipper) + bool preserve_intersection_clipper, + bool use_exact_kernel) { + if (LoopCGAL::verbose) { std::cout << "Cutting mesh with surface." << std::endl; } // Validate input meshes - if (!CGAL::is_valid_polygon_mesh(_mesh, LoopCGAL::verbose)) { + if (!CGAL::is_valid_polygon_mesh(_mesh, LoopCGAL::verbose)) + { std::cerr << "Error: Source mesh is invalid!" << std::endl; return; } - if (!CGAL::is_valid_polygon_mesh(clipper._mesh, LoopCGAL::verbose)) { + if (!CGAL::is_valid_polygon_mesh(clipper._mesh, LoopCGAL::verbose)) + { std::cerr << "Error: Clipper mesh is invalid!" << std::endl; return; } - if (_mesh.number_of_vertices() == 0 || _mesh.number_of_faces() == 0) { + if (_mesh.number_of_vertices() == 0 || _mesh.number_of_faces() == 0) + { std::cerr << "Error: Source mesh is empty!" << std::endl; return; } if (clipper._mesh.number_of_vertices() == 0 || - clipper._mesh.number_of_faces() == 0) { + clipper._mesh.number_of_faces() == 0) + { std::cerr << "Error: Clipper mesh is empty!" << std::endl; return; } @@ -321,24 +335,51 @@ void TriMesh::cutWithSurface(TriMesh &clipper, std::cout << "Clipping tm with clipper." << std::endl; } - try { - bool flag = - PMP::clip(_mesh, clipper._mesh, CGAL::parameters::clip_volume(false)); - - if (!flag) { + try + { + // bool flag = + // PMP::clip(_mesh, clipper._mesh, CGAL::parameters::clip_volume(false)); + bool flag = false; + try + { + if (use_exact_kernel){ + Exact_Mesh exact_clipper = convert_to_exact(clipper); + Exact_Mesh exact_mesh = convert_to_exact(*this); + flag = PMP::clip(exact_mesh, exact_clipper, CGAL::parameters::clip_volume(false)); + set_mesh(convert_to_double_mesh(exact_mesh)); + } + else{ + flag = PMP::clip(_mesh, clipper._mesh, CGAL::parameters::clip_volume(false)); + } + + } + catch (const std::exception &e) + { + std::cerr << "Corefinement failed: " << e.what() << std::endl; + } + if (!flag) + { std::cerr << "Warning: Clipping operation failed." << std::endl; - } else { - if (LoopCGAL::verbose) { + } + else + { + if (LoopCGAL::verbose) + { std::cout << "Clipping successful. Result has " << _mesh.number_of_vertices() << " vertices and " << _mesh.number_of_faces() << " faces." << std::endl; } } - } catch (const std::exception &e) { + } + catch (const std::exception &e) + { std::cerr << "Error during clipping: " << e.what() << std::endl; } - } else { - if (LoopCGAL::verbose) { + } + else + { + if (LoopCGAL::verbose) + { std::cout << "Meshes do not intersect. No clipping performed." << std::endl; } @@ -349,4 +390,316 @@ NumpyMesh TriMesh::save(double area_threshold, double duplicate_vertex_threshold) { return export_mesh(_mesh, area_threshold, duplicate_vertex_threshold); -} \ No newline at end of file +} + +void TriMesh::cut_with_implicit_function(const std::vector &property, double value, ImplicitCutMode cutmode) +{ + std::cout << "Cutting mesh with implicit function at value " << value << std::endl; + std::cout << "Mesh has " << _mesh.number_of_vertices() << " vertices and " + << _mesh.number_of_faces() << " faces." << std::endl; + std::cout << "Property size: " << property.size() << std::endl; + if (property.size() != _mesh.number_of_vertices()) + { + std::cerr << "Error: Property size does not match number of vertices." << std::endl; + return; + } + // Create a property map for vertex properties + typedef boost::property_map::type VertexIndexMap; + VertexIndexMap vim = get(boost::vertex_index, _mesh); + std::vector vertex_properties(_mesh.number_of_vertices()); + for (auto v : _mesh.vertices()) + { + vertex_properties[vim[v]] = property[vim[v]]; + } + auto property_map = boost::make_iterator_property_map( + vertex_properties.begin(), vim); + + // Build arrays similar to Python helper build_surface_arrays + // We'll need: edges map (ordered pair -> edge index), tri2edge mapping, edge->endpoints array + + // Map ordered vertex pair to an integer edge id + std::map, std::size_t> edge_index_map; + std::vector> edge_array; // endpoints + std::vector> tri_array; // triangles by vertex indices + + // Fill tri_array + for (auto f : _mesh.faces()) + { + std::array tri; + int i = 0; + for (auto v : vertices_around_face(_mesh.halfedge(f), _mesh)) + { + tri[i++] = vim[v]; + } + tri_array.push_back(tri); + } + + // Helper to get or create edge index + auto get_edge_id = [&](std::size_t a, std::size_t b) + { + if (a > b) + std::swap(a, b); + auto key = std::make_pair(a, b); + auto it = edge_index_map.find(key); + if (it != edge_index_map.end()) + return it->second; + std::size_t id = edge_array.size(); + edge_array.push_back(key); + edge_index_map[key] = id; + return id; + }; + + std::vector> tri2edge(tri_array.size()); + for (std::size_t ti = 0; ti < tri_array.size(); ++ti) + { + auto &tri = tri_array[ti]; + tri2edge[ti][0] = get_edge_id(tri[1], tri[2]); + tri2edge[ti][1] = get_edge_id(tri[2], tri[0]); + tri2edge[ti][2] = get_edge_id(tri[0], tri[1]); + } + + // Determine active triangles (all > value OR all < value OR any nan) + std::vector active(tri_array.size(), 0); + for (std::size_t ti = 0; ti < tri_array.size(); ++ti) + { + auto &tri = tri_array[ti]; + double v1 = vertex_properties[tri[0]]; + double v2 = vertex_properties[tri[1]]; + double v3 = vertex_properties[tri[2]]; + bool nan1 = std::isnan(v1); + bool nan2 = std::isnan(v2); + bool nan3 = std::isnan(v3); + if (nan1 || nan2 || nan3) + { + active[ti] = 1; + continue; + } + if ((v1 > value && v2 > value && v3 > value) || (v1 < value && v2 < value && v3 < value)) + { + active[ti] = 1; + } + } + + // Prepare new vertex list and new triangles similar to python + std::vector verts; + verts.reserve(_mesh.number_of_vertices()); + for (auto v : _mesh.vertices()) + verts.push_back(_mesh.point(v)); + std::vector newverts = verts; + std::vector newvals = vertex_properties; + + std::map new_point_on_edge; + std::vector> newtris(tri_array.begin(), tri_array.end()); + if (LoopCGAL::verbose) + { + std::cout << "Starting main loop over " << tri_array.size() << " triangles." << std::endl; + } + for (std::size_t t = 0; t < tri_array.size(); ++t) + { + if (active[t]) + { + continue; + } + auto tri = tri_array[t]; + // if all > value skip (hanging_wall in python) + if (vertex_properties[tri[0]] > value && vertex_properties[tri[1]] > value && vertex_properties[tri[2]] > value) + continue; + // for each edge of tri, check if edge crosses + for (auto eid : tri2edge[t]) + { + auto ends = edge_array[eid]; + double f0 = vertex_properties[ends.first]; + double f1 = vertex_properties[ends.second]; + if ((f0 > value && f1 > value) || (f0 < value && f1 < value)) + { + if (LoopCGAL::verbose) + { + std::cout << "Edge " << ends.first << "-" << ends.second << " does not cross the value " << value << std::endl; + } + continue; + } + double denom = (f1 - f0); + double ratio = 0.0; + if (std::abs(denom) < 1e-12) + ratio = 0.5; + else + ratio = (value - f0) / denom; + Point p0 = verts[ends.first]; + Point p1 = verts[ends.second]; + Point np = Point(p0.x() + ratio * (p1.x() - p0.x()), p0.y() + ratio * (p1.y() - p0.y()), p0.z() + ratio * (p1.z() - p0.z())); + newverts.push_back(np); + newvals.push_back(value); + new_point_on_edge[eid] = newverts.size() - 1; + } + + double v1 = vertex_properties[tri[0]]; + double v2 = vertex_properties[tri[1]]; + double v3 = vertex_properties[tri[2]]; + // replicate python cases + // convert tri to vector of 3 original indices and 2 new points + std::array extended = {tri[0], tri[1], tri[2], 0, 0}; + // retrieve relevant edges indices + std::size_t e01 = edge_index_map[std::make_pair(std::min(tri[0], tri[1]), std::max(tri[0], tri[1]))]; + std::size_t e12 = edge_index_map[std::make_pair(std::min(tri[1], tri[2]), std::max(tri[1], tri[2]))]; + std::size_t e20 = edge_index_map[std::make_pair(std::min(tri[2], tri[0]), std::max(tri[2], tri[0]))]; + // Get new points where available + std::size_t np_e01 = new_point_on_edge.count(e01) ? new_point_on_edge[e01] : SIZE_MAX; + std::size_t np_e12 = new_point_on_edge.count(e12) ? new_point_on_edge[e12] : SIZE_MAX; + std::size_t np_e20 = new_point_on_edge.count(e20) ? new_point_on_edge[e20] : SIZE_MAX; + // Helper to append triangle + auto append_tri = [&](std::array tarr) + { newtris.push_back(tarr); }; + + // CASE 1: v1 > value and v2 > value and v3 value && v2 > value && v3 < value) + { + std::size_t p1 = np_e12; + std::size_t p2 = np_e20; + extended[3] = p1; + extended[4] = p2; + std::array m1 = {extended[0], extended[1], extended[3]}; + std::array m2 = {extended[0], extended[3], extended[4]}; + std::array m3 = {extended[4], extended[3], extended[2]}; + newtris[t] = m1; + append_tri(m2); + append_tri(m3); + if (LoopCGAL::verbose) + { + std::cout << "CASE 1 executed" << std::endl; + } + continue; + } + // CASE 2 + if (v1 > value && v2 < value && v3 > value) + { + std::size_t p1 = np_e01; + std::size_t p2 = np_e12; + extended[3] = p1; + extended[4] = p2; + std::array m1 = {extended[0], extended[3], extended[2]}; + std::array m2 = {extended[3], extended[4], extended[2]}; + std::array m3 = {extended[3], extended[1], extended[4]}; + newtris[t] = m1; + append_tri(m2); + append_tri(m3); + if (LoopCGAL::verbose) + { + std::cout << "CASE 2 executed" << std::endl; + } + continue; + } + // CASE 3 + if (v1 < value && v2 > value && v3 > value) + { + std::size_t p1 = np_e01; + std::size_t p2 = np_e20; + extended[3] = p1; + extended[4] = p2; + std::array m1 = {extended[0], extended[3], extended[4]}; + std::array m2 = {extended[3], extended[1], extended[2]}; + std::array m3 = {extended[4], extended[3], extended[2]}; + newtris[t] = m1; + append_tri(m2); + append_tri(m3); + if (LoopCGAL::verbose) + { + std::cout << "CASE 3 executed" << std::endl; + } + continue; + } + // CASE 5 + if (v1 < value && v2 < value && v3 > value) + { + std::size_t p1 = np_e12; + std::size_t p2 = np_e20; + extended[3] = p1; + extended[4] = p2; + std::array m1 = {extended[0], extended[1], extended[3]}; + std::array m2 = {extended[0], extended[3], extended[4]}; + std::array m3 = {extended[4], extended[3], extended[2]}; + newtris[t] = m1; + append_tri(m2); + append_tri(m3); + if (LoopCGAL::verbose) + { + std::cout << "CASE 5 executed" << std::endl; + } + continue; + } + // CASE 6 + if (v1 < value && v2 > value && v3 < value) + { + std::size_t p1 = np_e01; + std::size_t p2 = np_e12; + extended[3] = p1; + extended[4] = p2; + std::array m1 = {extended[0], extended[3], extended[2]}; + std::array m2 = {extended[3], extended[4], extended[2]}; + std::array m3 = {extended[3], extended[1], extended[4]}; + newtris[t] = m1; + append_tri(m2); + append_tri(m3); + if (LoopCGAL::verbose) + { + std::cout << "CASE 6 executed" << std::endl; + } + continue; + } + // CASE 7 + if (v1 > value && v2 < value && v3 < value) + { + std::size_t p1 = np_e01; + std::size_t p2 = np_e20; + extended[3] = p1; + extended[4] = p2; + std::array m1 = {extended[0], extended[3], extended[4]}; + std::array m2 = {extended[3], extended[2], extended[4]}; + std::array m3 = {extended[3], extended[1], extended[2]}; + newtris[t] = m1; + append_tri(m2); + append_tri(m3); + if (LoopCGAL::verbose) + { + std::cout << "CASE 7 executed" << std::endl; + } + continue; + } + } + + // Build new CGAL mesh from newverts and newtris + TriangleMesh newmesh; + std::vector new_vhandles; + new_vhandles.reserve(newverts.size()); + for (auto &p : newverts) + new_vhandles.push_back(newmesh.add_vertex(p)); + for (auto &tri : newtris) + { + // skip degenerate + if (tri[0] == tri[1] || tri[1] == tri[2] || tri[0] == tri[2]) + continue; + if (ImplicitCutMode::KEEP_NEGATIVE_SIDE == cutmode) + { + double v0 = newvals[tri[0]]; + double v1 = newvals[tri[1]]; + double v2 = newvals[tri[2]]; + if (v0 > value && v1 > value && v2 > value) + { + continue; + } + } + if (ImplicitCutMode::KEEP_POSITIVE_SIDE == cutmode) + { + double v0 = newvals[tri[0]]; + double v1 = newvals[tri[1]]; + double v2 = newvals[tri[2]]; + if (v0 < value && v1 < value && v2 < value) + { + continue; + } + } + newmesh.add_face(new_vhandles[tri[0]], new_vhandles[tri[1]], new_vhandles[tri[2]]); + } + + // Replace internal mesh + _mesh = std::move(newmesh); +} diff --git a/src/mesh.h b/src/mesh.h index ca4c3fc..c491b7e 100644 --- a/src/mesh.h +++ b/src/mesh.h @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -10,6 +11,9 @@ #include #include // For std::pair #include +#include "meshenums.h" +typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_K; +typedef CGAL::Surface_mesh Exact_Mesh; typedef CGAL::Simple_cartesian Kernel; typedef Kernel::Point_3 Point; typedef CGAL::Surface_mesh TriangleMesh; @@ -27,18 +31,21 @@ class TriMesh // Method to cut the mesh with another surface object void cutWithSurface(TriMesh &surface, bool preserve_intersection = false, - bool preserve_intersection_clipper = false); + bool preserve_intersection_clipper = false, + bool use_exact_kernel = true); // Method to remesh the triangle mesh void remesh(bool split_long_edges, double target_edge_length, int number_of_iterations, bool protect_constraints, bool relax_constraints); void init(); + void cut_with_implicit_function(const std::vector& property, double value, ImplicitCutMode cutmode = ImplicitCutMode::KEEP_POSITIVE_SIDE); // Getters for mesh properties void reverseFaceOrientation(); NumpyMesh save(double area_threshold, double duplicate_vertex_threshold); void add_fixed_edges(const pybind11::array_t &pairs); - + const TriangleMesh& get_mesh() const { return _mesh; } + void set_mesh(const TriangleMesh& mesh) { _mesh = mesh; } private: std::set _fixedEdges; TriangleMesh _mesh; // The underlying CGAL surface mesh diff --git a/src/meshenums.h b/src/meshenums.h new file mode 100644 index 0000000..357bbc5 --- /dev/null +++ b/src/meshenums.h @@ -0,0 +1,9 @@ +#ifndef MESH_ENUMS_H +#define MESH_ENUMS_H +enum class ImplicitCutMode +{ + PRESERVE_INTERSECTION, + KEEP_POSITIVE_SIDE, + KEEP_NEGATIVE_SIDE +}; +#endif // MESH_ENUMS_H \ No newline at end of file diff --git a/src/meshutils.cpp b/src/meshutils.cpp index 63e2500..7bffc50 100644 --- a/src/meshutils.cpp +++ b/src/meshutils.cpp @@ -256,5 +256,45 @@ NumpyMesh export_mesh(const TriangleMesh &tm, double area_threshold, NumpyMesh result; result.vertices = vertices_array; result.triangles = triangles_array; + return result; +} +Exact_Mesh convert_to_exact(const TriMesh& input) { + Exact_Mesh result; + std::map vmap; + + for (auto v : vertices(input.get_mesh())) { + const auto& p = input.get_mesh().point(v); + Exact_K::Point_3 ep(p.x(), p.y(), p.z()); + vmap[v] = result.add_vertex(ep); + } + + for (auto f : faces(input.get_mesh())) { + std::vector face_vertices; + for (auto v : vertices_around_face(input.get_mesh().halfedge(f), input.get_mesh())) { + face_vertices.push_back(vmap[v]); + } + result.add_face(face_vertices); + } + + return result; +} +TriangleMesh convert_to_double_mesh(const Exact_Mesh& input) { + TriangleMesh result; + std::map vmap; + + for (auto v : vertices(input)) { + const auto& p = input.point(v); + Point dp(CGAL::to_double(p.x()), CGAL::to_double(p.y()), CGAL::to_double(p.z())); + vmap[v] = result.add_vertex(dp); + } + + for (auto f : faces(input)) { + std::vector face_vertices; + for (auto v : vertices_around_face(input.halfedge(f), input)) { + face_vertices.push_back(vmap[v]); + } + result.add_face(face_vertices); + } + return result; } \ No newline at end of file diff --git a/src/meshutils.h b/src/meshutils.h index c4fc9aa..5307f45 100644 --- a/src/meshutils.h +++ b/src/meshutils.h @@ -8,4 +8,7 @@ NumpyMesh export_mesh(const TriangleMesh &tm, double area_threshold, double calculate_triangle_area(const std::array &v1, const std::array &v2, const std::array &v3); +Exact_Mesh convert_to_exact(const TriMesh& input); +TriangleMesh convert_to_double_mesh(const Exact_Mesh& input); + #endif // MESHUTILS_H diff --git a/tests/test_degenerate.py b/tests/test_degenerate.py index 0bd60e4..f263c5d 100644 --- a/tests/test_degenerate.py +++ b/tests/test_degenerate.py @@ -1,22 +1,30 @@ import numpy as np -import loop_cgal import pyvista as pv -def test_degenerate_triangles(): - """Test handling of degenerate triangles in TriMesh.""" - # Create a surface with degenerate triangles - points = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0.5, 0.5, 0]]) - faces = np.array([[0, 1, 2], [2, 3, 4], [4, 0, 0]]) # Degenerate triangle (4 is not a valid vertex) - - surface = pv.PolyData.from_regular_faces(points, faces) - - try: - tri_mesh = loop_cgal.TriMesh(surface) - print("TriMesh created successfully with degenerate triangles.") - except ValueError as e: - print(f"ValueError: {e}") - - -if __name__ == "__main__": - loop_cgal.set_verbose(True) - test_degenerate_triangles() - print("Test completed.") \ No newline at end of file +import pytest +import loop_cgal + + +def make_invalid_polydata_invalid_index(): + # face references a non-existent vertex index (4) + points = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]]) + faces = np.array([[0, 1, 2], [2, 3, 4]]) + return pv.PolyData.from_regular_faces(points, faces) + + +def make_polydata_degenerate_triangle(): + # triangle repeats a vertex index -> degenerate + points = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0]]) + faces = np.array([[0, 1, 1]]) + return pv.PolyData.from_regular_faces(points, faces) + + +def test_invalid_index_raises_value_error(): + surface = make_invalid_polydata_invalid_index() + with pytest.raises(ValueError, match="exceed vertex count"): + _ = loop_cgal.TriMesh(surface) + + +def test_degenerate_triangle_raises_value_error(): + surface = make_polydata_degenerate_triangle() + with pytest.raises(ValueError, match="degenerate triangles"): + _ = loop_cgal.TriMesh(surface) \ No newline at end of file diff --git a/tests/test_mesh_operations.py b/tests/test_mesh_operations.py new file mode 100644 index 0000000..4eb9c5d --- /dev/null +++ b/tests/test_mesh_operations.py @@ -0,0 +1,75 @@ +from __future__ import annotations +import numpy as np +import pyvista as pv +import pytest +import loop_cgal +from loop_cgal._loop_cgal import ImplicitCutMode + + +@pytest.fixture +def square_surface(): + # Unit square made of two triangles + return pv.Plane(center=(0,0,0),direction=(0,0,1),i_size=1.0,j_size=1.0) + + +@pytest.fixture +def clipper_surface(): + # A square that overlaps half of the unit square + return pv.Plane(center=(0,0,0),direction=(1,0,0),i_size=2.0,j_size=2.0) + + +# @pytest.mark.parametrize("remesh_kwargs", [ +# {"split_long_edges": True, "target_edge_length": 0.2, "number_of_iterations": 1, "protect_constraints": True, "relax_constraints": False}, +# {"split_long_edges": False, "target_edge_length": 0.02, "number_of_iterations": 2, "protect_constraints": True, "relax_constraints": True}, +# ]) +def test_loading_and_saving(square_surface): + tm = loop_cgal.TriMesh(square_surface) + saved = tm.save() + verts = np.array(saved.vertices) + tris = np.array(saved.triangles) + assert verts.ndim == 2 and verts.shape[1] == 3 + assert tris.ndim == 2 and tris.shape[1] == 3 + assert verts.shape[0] > 0 + assert tris.shape[0] > 0 + + +def test_cut_with_surface(square_surface, clipper_surface): + tm = loop_cgal.TriMesh(square_surface) + clip = loop_cgal.TriMesh(clipper_surface) + before = np.array(tm.save().triangles).shape[0] + tm.cut_with_surface(clip) + after = np.array(tm.save().triangles).shape[0] + # If clipper intersects, faces should be non-zero and not increase + assert after >= 0 + assert after <= before + + +@pytest.mark.parametrize("kwargs", [ + {"split_long_edges": True, "target_edge_length": 0.25, "number_of_iterations": 1, "protect_constraints": True, "relax_constraints": False}, + {"split_long_edges": True, "target_edge_length": 0.05, "number_of_iterations": 2, "protect_constraints": False, "relax_constraints": True}, +]) +def test_remesh_changes_vertices(square_surface, kwargs): + tm = loop_cgal.TriMesh(square_surface) + + # Call remesh using keyword args compatible with the binding + tm.remesh(kwargs["split_long_edges"], kwargs["target_edge_length"], kwargs["number_of_iterations"], kwargs["protect_constraints"], kwargs["relax_constraints"]) + after_v = np.array(tm.save().vertices).shape[0] + # Remesh should produce a valid mesh + assert after_v > 0 + # Either vertices increase due to splitting or stay similar; ensure no catastrophic collapse + # assert after_v >= 0.5 * before_v + + +def test_cut_with_implicit_function(square_surface): + tm = loop_cgal.TriMesh(square_surface) + # create a scalar property that varies across vertices + saved = tm.save() + nverts = np.array(saved.vertices).shape[0] + prop = [float(i) / max(1, (nverts - 1)) for i in range(nverts)] + # cut at 0.5 keeping positive side + tm.cut_with_implicit_function(prop, 0.5, ImplicitCutMode.KEEP_POSITIVE_SIDE) + res = tm.save() + v = np.array(res.vertices).shape[0] + f = np.array(res.triangles).shape[0] + assert v >= 0 + assert f >= 0