|
5 | 5 | from __future__ import annotations |
6 | 6 |
|
7 | 7 | import logging |
8 | | -from typing import cast |
| 8 | +import os |
| 9 | +from typing import Optional, cast |
9 | 10 |
|
10 | 11 | import itk |
11 | 12 | import numpy as np |
@@ -70,6 +71,84 @@ def extract_contours( |
70 | 71 |
|
71 | 72 | return contours |
72 | 73 |
|
| 74 | + def extract_mesh( |
| 75 | + self, surface: pv.PolyData, mesh_target_reduction: float = 0.0 |
| 76 | + ) -> Optional[pv.UnstructuredGrid]: |
| 77 | + """Generate a tetrahedral volume mesh (VTU) from a closed surface via netgen. |
| 78 | +
|
| 79 | + Optionally decimates the surface with |
| 80 | + :meth:`pyvista.PolyDataFilters.decimate_pro` before meshing — netgen |
| 81 | + has no post-hoc decimation of its own, so a coarser input surface is |
| 82 | + the only way to get a coarser tetrahedral mesh out. |
| 83 | +
|
| 84 | + Builds netgen's surface mesh directly from *surface*'s indexed |
| 85 | + points/triangles rather than round-tripping through an STL file. |
| 86 | + STL has no shared-vertex topology, so every triangle corner is |
| 87 | + written independently; float32 rounding during that round-trip can |
| 88 | + make two corners that were exactly the same point diverge by ~1e-6, |
| 89 | + which is well within netgen's own "identical point" merge tolerance |
| 90 | + and makes its STL reader spin forever trying to reconcile them |
| 91 | + (observed hang on Taubin-smoothed surfaces, e.g. from |
| 92 | + :meth:`extract_contours`). Adding points/triangles directly preserves |
| 93 | + the exact shared-vertex indices already present in *surface*, so |
| 94 | + there is nothing for netgen to reconcile. |
| 95 | +
|
| 96 | + Args: |
| 97 | + surface: Closed, triangulated surface for one anatomy group (as |
| 98 | + returned by :meth:`extract_contours`). |
| 99 | + mesh_target_reduction: Fraction in ``[0, 1)`` of surface triangles |
| 100 | + to remove via ``decimate_pro(mesh_target_reduction, |
| 101 | + preserve_topology=True)`` before meshing. ``0.0`` (default) |
| 102 | + skips decimation and meshes the surface as given. |
| 103 | +
|
| 104 | + Returns: |
| 105 | + :class:`pyvista.UnstructuredGrid` of tetrahedral cells, or |
| 106 | + ``None`` if the surface is empty or netgen produced no volume |
| 107 | + mesh from it. |
| 108 | + """ |
| 109 | + if surface.n_points == 0: |
| 110 | + return None |
| 111 | + |
| 112 | + meshing_surface = surface.triangulate().clean() |
| 113 | + if mesh_target_reduction > 0.0: |
| 114 | + meshing_surface = meshing_surface.decimate_pro( |
| 115 | + mesh_target_reduction, preserve_topology=True |
| 116 | + ) |
| 117 | + |
| 118 | + import netgen.meshing as ngm # noqa: PLC0415 |
| 119 | + |
| 120 | + triangles = meshing_surface.faces.reshape(-1, 4)[:, 1:4] |
| 121 | + ngmesh = ngm.Mesh() |
| 122 | + ngmesh.dim = 3 |
| 123 | + point_ids = [ |
| 124 | + ngmesh.Add(ngm.MeshPoint(ngm.Pnt(*point))) |
| 125 | + for point in meshing_surface.points |
| 126 | + ] |
| 127 | + face_descriptor = ngmesh.Add(ngm.FaceDescriptor(surfnr=1, domin=1, domout=0)) |
| 128 | + for triangle in triangles: |
| 129 | + ngmesh.Add(ngm.Element2D(face_descriptor, [point_ids[i] for i in triangle])) |
| 130 | + ngmesh.GenerateVolumeMesh() |
| 131 | + |
| 132 | + elements = ngmesh.Elements3D() |
| 133 | + if len(elements) == 0: |
| 134 | + self.log_warning("netgen produced no volume mesh for this surface") |
| 135 | + return None |
| 136 | + |
| 137 | + points = ngmesh.Coordinates() |
| 138 | + tets = np.array( |
| 139 | + [[vertex.nr - 1 for vertex in element.vertices] for element in elements], |
| 140 | + dtype=np.int64, |
| 141 | + ) |
| 142 | + # netgen's tet vertex order is opposite VTK_TETRA's right-hand |
| 143 | + # convention (confirmed by negative cell volumes); swapping the last |
| 144 | + # two indices restores positive-volume orientation. |
| 145 | + tets = tets[:, [0, 1, 3, 2]] |
| 146 | + cells = np.hstack( |
| 147 | + [np.full((tets.shape[0], 1), 4, dtype=np.int64), tets] |
| 148 | + ).flatten() |
| 149 | + cell_types = np.full(tets.shape[0], pv.CellType.TETRA, dtype=np.uint8) |
| 150 | + return pv.UnstructuredGrid(cells, cell_types, points) |
| 151 | + |
73 | 152 | def transform_contours( |
74 | 153 | self, |
75 | 154 | contours: pv.PolyData, |
@@ -462,3 +541,130 @@ def create_deformation_field( |
462 | 541 | ) |
463 | 542 |
|
464 | 543 | return deformation_field_img |
| 544 | + |
| 545 | + # ─────────────────────────── I/O helpers ─────────────────────────────── |
| 546 | + |
| 547 | + @staticmethod |
| 548 | + def save_surfaces( |
| 549 | + surfaces: dict[str, pv.PolyData], |
| 550 | + output_dir: str, |
| 551 | + prefix: str = "", |
| 552 | + ) -> dict[str, str]: |
| 553 | + """Save each named surface to its own VTP file. |
| 554 | +
|
| 555 | + Args: |
| 556 | + surfaces: Mapping of name → surface (e.g. the ``'surfaces'`` |
| 557 | + value from :meth:`WorkflowConvertImageToVTK.process`). |
| 558 | + output_dir: Directory to write files into (created if absent). |
| 559 | + prefix: Optional filename prefix. Each file is named |
| 560 | + ``{prefix}_{name}.vtp`` (or ``{name}.vtp`` when *prefix* is empty). |
| 561 | +
|
| 562 | + Returns: |
| 563 | + Mapping of name → absolute path of the saved file. |
| 564 | + """ |
| 565 | + os.makedirs(output_dir, exist_ok=True) |
| 566 | + saved: dict[str, str] = {} |
| 567 | + for name, surface in surfaces.items(): |
| 568 | + stem = f"{prefix}_{name}" if prefix else name |
| 569 | + path = os.path.join(output_dir, f"{stem}.vtp") |
| 570 | + surface.save(path) |
| 571 | + saved[name] = path |
| 572 | + return saved |
| 573 | + |
| 574 | + @staticmethod |
| 575 | + def save_meshes( |
| 576 | + meshes: dict[str, pv.UnstructuredGrid], |
| 577 | + output_dir: str, |
| 578 | + prefix: str = "", |
| 579 | + ) -> dict[str, str]: |
| 580 | + """Save each named volume mesh to its own VTU file. |
| 581 | +
|
| 582 | + Args: |
| 583 | + meshes: Mapping of name → mesh (e.g. the ``'meshes'`` value from |
| 584 | + :meth:`WorkflowConvertImageToVTK.process`). |
| 585 | + output_dir: Directory to write files into (created if absent). |
| 586 | + prefix: Optional filename prefix. Each file is named |
| 587 | + ``{prefix}_{name}.vtu`` (or ``{name}.vtu`` when *prefix* is empty). |
| 588 | +
|
| 589 | + Returns: |
| 590 | + Mapping of name → absolute path of the saved file. |
| 591 | + """ |
| 592 | + os.makedirs(output_dir, exist_ok=True) |
| 593 | + saved: dict[str, str] = {} |
| 594 | + for name, mesh in meshes.items(): |
| 595 | + stem = f"{prefix}_{name}" if prefix else name |
| 596 | + path = os.path.join(output_dir, f"{stem}.vtu") |
| 597 | + mesh.save(path) |
| 598 | + saved[name] = path |
| 599 | + return saved |
| 600 | + |
| 601 | + @staticmethod |
| 602 | + def save_combined_surface( |
| 603 | + surfaces: dict[str, pv.PolyData], |
| 604 | + output_dir: str, |
| 605 | + prefix: str = "", |
| 606 | + ) -> str: |
| 607 | + """Merge all named surfaces into a single VTP file. |
| 608 | +
|
| 609 | + The merged mesh retains per-cell ``Color`` (RGBA uint8) from each |
| 610 | + surface's annotation, enabling colour-by-anatomy rendering in |
| 611 | + Paraview, PyVista, etc. Per-object ``field_data`` is not preserved |
| 612 | + in the merged file. |
| 613 | +
|
| 614 | + Args: |
| 615 | + surfaces: Mapping of name → surface. |
| 616 | + output_dir: Directory to write the file into (created if absent). |
| 617 | + prefix: Optional filename prefix. Output is ``{prefix}_surfaces.vtp`` |
| 618 | + (or ``surfaces.vtp`` when *prefix* is empty). |
| 619 | +
|
| 620 | + Returns: |
| 621 | + Absolute path to the saved VTP file. |
| 622 | +
|
| 623 | + Raises: |
| 624 | + ValueError: If *surfaces* is empty. |
| 625 | + """ |
| 626 | + if not surfaces: |
| 627 | + raise ValueError("No surfaces to save.") |
| 628 | + os.makedirs(output_dir, exist_ok=True) |
| 629 | + stem = f"{prefix}_surfaces" if prefix else "surfaces" |
| 630 | + output_file = os.path.join(output_dir, f"{stem}.vtp") |
| 631 | + merged = cast( |
| 632 | + pv.PolyData, pv.merge(list(surfaces.values()), merge_points=False) |
| 633 | + ) |
| 634 | + merged.save(output_file) |
| 635 | + return output_file |
| 636 | + |
| 637 | + @staticmethod |
| 638 | + def save_combined_mesh( |
| 639 | + meshes: dict[str, pv.UnstructuredGrid], |
| 640 | + output_dir: str, |
| 641 | + prefix: str = "", |
| 642 | + ) -> str: |
| 643 | + """Merge all named volume meshes into a single VTU file. |
| 644 | +
|
| 645 | + The merged mesh retains per-cell ``Color`` (RGBA uint8) from each |
| 646 | + mesh's annotation. Per-object ``field_data`` is not preserved in the |
| 647 | + merged file. |
| 648 | +
|
| 649 | + Args: |
| 650 | + meshes: Mapping of name → volume mesh. |
| 651 | + output_dir: Directory to write the file into (created if absent). |
| 652 | + prefix: Optional filename prefix. Output is ``{prefix}_meshes.vtu`` |
| 653 | + (or ``meshes.vtu`` when *prefix* is empty). |
| 654 | +
|
| 655 | + Returns: |
| 656 | + Absolute path to the saved VTU file. |
| 657 | +
|
| 658 | + Raises: |
| 659 | + ValueError: If *meshes* is empty. |
| 660 | + """ |
| 661 | + if not meshes: |
| 662 | + raise ValueError("No meshes to save.") |
| 663 | + os.makedirs(output_dir, exist_ok=True) |
| 664 | + stem = f"{prefix}_meshes" if prefix else "meshes" |
| 665 | + output_file = os.path.join(output_dir, f"{stem}.vtu") |
| 666 | + merged = cast( |
| 667 | + pv.UnstructuredGrid, pv.merge(list(meshes.values()), merge_points=False) |
| 668 | + ) |
| 669 | + merged.save(output_file) |
| 670 | + return output_file |
0 commit comments