Volume remeshing

This code implements a variation of the volume meshing algorithm described in "Convex Polyhedral Meshing for Robust Solid Modeling" by Lorenzo Diazzi and Marco Attene (ACM Trans Graphics Vol 40, N. 6, Procs of SIGGRAPH Asia 2021). You may download a copy of the paper preprint here: http://arxiv.org/abs/2109.14434 The original code is available here: https://github.com/MarcoAttene/VolumeMesher

This fork adds, on top of the original algorithm:

• Optional, combinable inputs. A run takes any non-empty combination of an input surface (triangles), a set of edges, and a set of points — see Inserting edges and points.
• Exact surface tracking. For each output element you can recover which input primitive it came from: output faces → input triangles, output edges → input edges, output vertices → input points — see Tracking output back to the input.
• Embedding a surface (and edges/points) in an existing tetrahedral mesh, programmatically via embed.h.
• A 2D mode (-2d): the arrangement of a soup of 2D segments, output as a triangulation with the same kind of exact provenance — see 2D mode.
• Python bindings (pip install volume-remesher): every option of the command line tool as a function, with numpy arrays in and out and files that are byte-identical to the tool's — see Python.

Building

Clone this repository and configure with CMake:

git clone https://github.com/MarcoAttene/VolumeRemesher
cd VolumeRemesher
mkdir build
cd build
cmake ..

This produces a build configuration for your system (a mesh_generator.sln on Windows/MSVC, a Makefile on Linux/macOS). Build the mesh_generator target as usual. The integration tests are built by default (VOLUMEREMESHER_BUILD_TESTS=ON) — run them with ctest, or configure with -DVOLUMEREMESHER_BUILD_TESTS=OFF to skip them. (Some heavy tracking tests only run in a Release build.)

The legacy graph-cut cell classifier (-c graphcut, see Inside/outside classification) is not built by default because its third-party code is not LGPL compatible; -DVOLUMEREMESHER_WITH_GRAPHCUT=ON compiles it in, and ctest then also runs the suite against it.

Python

Installation

Wheels for Linux (x86-64, aarch64), macOS (Apple Silicon) and Windows (x64), for CPython 3.9 and later:

pip install volume-remesher

The x86-64 wheels are compiled with AVX2/FMA, like the command line tool: they need a CPU from 2013 or later (Intel Haswell, AMD Excavator). On an older CPU, or on another platform, build from source instead; this needs a C++20 compiler and CMake 3.25 (both are what the tool itself needs), and compiles the same code with the same flags:

pip install git+https://github.com/wildmeshing/VolumeRemesher   # or, from a checkout: pip install .

To work on the package in-tree, configure CMake with -DVOLUMEREMESHER_BUILD_PYTHON=ON (after pip install nanobind numpy): the module is placed in python/volume_remesher/, and PYTHONPATH=python python -m pytest tests/python runs the Python test suite against the mesh_generator of the same build.

Usage

Every input is a path or numpy arrays; every result is a frozen dataclass of numpy arrays with a save() that writes exactly the bytes mesh_generator writes (the test suite checks this on every model and option):

import volume_remesher as vr

# 1) Volume mesh (polyhedral) of a closed surface      mesh_generator models/bunny.off
mesh = vr.remesh("models/bunny.off")
mesh.vertices            # (n, 3) float64
mesh.cells               # the polyhedral cells, each a list of vertex indices
mesh.save("volume.msh")

# 2) Tetrahedral mesh of the enclosed solid            mesh_generator models/bunny.off -t
mesh = vr.remesh("models/bunny.off", tetrahedralize=True)
mesh.tets                # (m, 4) uint32
mesh.save("volume.tet")  # + volume.tet.triangleprov, volume.tet.groups

# 3) Boolean union of two overlapping surfaces         mesh_generator cube_subdiv.off U sphere.off
mesh = vr.boolean("models/cube_subdiv.off", "models/sphere.off", "union")   # or "intersection", "difference"

# 4) Exact surface tracking of a tetrahedralization    mesh_generator cube_subdiv.off -t -r -a
mesh = vr.remesh("models/cube_subdiv.off", tetrahedralize=True, exact=True, keep_all_cells=True)
mesh.exact_vertices      # list of (Fraction, Fraction, Fraction): the exact positions
mesh.triangle_provenance # per coplanar group of input triangles, the {tet, v0, v1, v2} faces tiling it
mesh.constraint_groups   # (c, 2): input triangle -> coplanar group
# Cheaper variants, Python only: face_provenance="faces" gives just which tet faces lie on the
# input (mesh.surface_faces, a per-tet bitmask) without the coplanar grouping, "none" nothing;
# round_vertices=True makes mesh.vertices the nearest doubles of the exact positions.

# 5) Insert edges and points, with tracking            mesh_generator sphere.off -t -r -a -e ... -p ...
mesh = vr.remesh("models/sphere.off", edges="models/sphere_edges.off", points="models/sphere_points.off",
                 tetrahedralize=True, exact=True, keep_all_cells=True)
mesh.edge_provenance[i]  # the {tet, v0, v1} output edges tiling input edge mesh.edge_ids[i]
mesh.point_provenance    # (p, 3): {point_id, tet, output vertex}, -1 where the point did not survive

# 6) Embed a surface into an existing tetrahedral mesh  mesh_generator surface.off U mesh.tet
mesh = vr.embed("models/Octocat.off", "models/Octocat.bg.tet", tetrahedralize=True)

# 7) 2D: arrangement of a segment soup                 mesh_generator models/2d/england.obj -2d
arr = vr.arrange_2d("models/2d/england.obj")
arr.vertices, arr.triangles, arr.segment_provenance
arr.save("triangulation.off")   # + .segmentprov, .pointprov

Working with arrays

Arrays work wherever a path does: vr.remesh((vertices, triangles), ...), edges=(edge_vertices, edge_pairs), points=point_array, vr.embed(surface, (tet_vertices, tets)), vr.arrange_2d((points_2d, segment_pairs)). Coordinates are float64 (n, 3) (or (n, 2) in 2D) and indices any integer array, (m, 3) for triangles, (m, 4) for tets, (m, 2) for edges and segments. The readers are available on their own (vr.read_off, vr.read_edges_off, vr.read_points_off, vr.read_tet_mesh, vr.read_obj_segments) and parse exactly as the tool does.

The arrangement of two crossing segments, the diagonals of the square [0,2]², from arrays:

import numpy as np
import volume_remesher as vr

points = np.array([[0.0, 0.0], [2.0, 2.0], [0.0, 2.0], [2.0, 0.0]])
segments = np.array([[0, 1], [2, 3]])          # (0,0)-(2,2) and (0,2)-(2,0)

arr = vr.arrange_2d((points, segments), exact=True)
arr                       # Arrangement(9 vertices, 12 triangles)
arr.vertices              # (9, 2): the 4 input points, the crossing at (1, 1), and the
                          #         4 corners of the bounding box expanded by 10%
arr.triangles             # (12, 3) uint32, counterclockwise
arr.exact_vertices[1]     # (Fraction(1, 1), Fraction(1, 1)): the crossing, exactly

# Each input segment is exactly the union of output edges, in order from its first endpoint:
arr.segment_provenance[0] # [[2 4 1]      triangle 2 has edge (4, 1): (0,0) -> (1,1)
                          #  [0 1 2]]     triangle 0 has edge (1, 2): (1,1) -> (2,2)
arr.point_provenance      # (4, 2): for each input point, a triangle containing it and its
                          #         output vertex index
arr.save("x.off")         # the same OFF (+ .segmentprov, .pointprov, .rational) as the tool

And a 3D example, the tetrahedral mesh of a unit cube given as arrays, then its union with a shifted copy:

V = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0],
              [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]], dtype=float)
F = np.array([[0, 2, 1], [0, 3, 2], [4, 5, 6], [4, 6, 7], [0, 1, 5], [0, 5, 4],
              [2, 3, 7], [2, 7, 6], [1, 2, 6], [1, 6, 5], [0, 4, 7], [0, 7, 3]])

mesh = vr.remesh((V, F), tetrahedralize=True)          # VolumeMesh(8 vertices, 6 tets)
both = vr.boolean((V, F), (V + 0.5, F), "union", tetrahedralize=True)   # VolumeMesh(22 vertices, 42 tets)
both.vertices, both.tets                                # numpy arrays, ready for your own pipeline

Errors — a missing or malformed file, a vertex outside the tet mesh, an input that cannot be meshed — raise vr.VolumeRemesherError with the tool's message; the interpreter is never terminated. The pipeline runs on its own thread with the GIL released, and its floating point rounding mode never leaks into the calling thread. Its parallel passes use vr.get_num_threads() threads: vr.set_num_threads(n), else the VOLUMEREMESHER_NUM_THREADS environment variable, else the CPUs the process may run on (its affinity mask on Linux, e.g. under taskset or a batch scheduler; every CPU of the machine on macOS and Windows). Progress output (verbose=True) goes to the C standard output, as in the tool; log=True appends to mesh_generator.log in the working directory.

The command line, in Python

python -m volume_remesher (also installed as volume-remesher) takes mesh_generator's exact arguments and writes its exact files, with the same exit status and messages:

python -m volume_remesher models/bunny.off -t -r      # == mesh_generator models/bunny.off -t -r

Performance is that of the tool: the same compiled code runs, and the Python side only copies the result into its arrays (on the bundled models, tests/python/bench.py measures the Python call at 0.95-1.0x the tool's wall time, slightly less without the file write).

Command-line usage

mesh_generator [options] [inputfile_A.off] [bool_opcode inputfile_B.off] [-e edges.off] [-p points.off]

The three inputs — the surface inputfile_A.off, the -e edges and the -p points — are all optional; provide any non-empty combination (only a surface, only edges, only points, or any mix). Run mesh_generator with no arguments to print the built-in help.

Options

Inside/outside classification

Unless -a is given, the cells of the arrangement are classified as inside or outside the input surface and only the inside ones are kept (per surface for a boolean, then combined by the opcode). The default classifier, -c winding, is the winding number of the surface:

• For a closed, consistently oriented input the winding number is computed exactly on the arrangement itself: it is constant on every cell, equal across a face that no input triangle covers and changing by the signed number of covering triangles across one that some do, so a walk from the outside assigns it to every cell with exact orientation predicates only. A cell is inside iff its winding number is nonzero (so a self-intersecting input is the union of what its shells enclose). Before the walk the orientation of the triangles is made consistent across manifold edges and each closed component is oriented outward, which changes nothing for a well-formed input and makes the result meaningful for one with flipped faces.
• For an input that is open or that cannot be oriented consistently (the walk detects this) every cell takes the generalized winding number of the surface at a point inside it (Jacobson et al. 2013, evaluated with the hierarchical expansion of Barill et al. 2018), and is inside iff that exceeds 1/2 in magnitude, as fTetWild filters its tetrahedra. The evaluation uses only IEEE arithmetic (with a deterministic atan2), so the output stays byte-identical across platforms. Capping the input first (-f) keeps it on the exact path.

-c graphcut is the classifier the tool used before: an area-weighted min-cut over the dual graph of the cells. It gives the same result as the winding number on closed, consistently oriented inputs (the test suite checks this cell by cell) and a different, heuristic closure on open or inconsistently oriented ones. It depends on the GCO and maxflow libraries in third_party/graph_cut, whose licenses (research use only; GPL) are not compatible with this library's LGPL, so it is compiled only on request (-DVOLUMEREMESHER_WITH_GRAPHCUT=ON) and never shipped in the Python wheels.

Input formats

• Surface (inputfile_A.off, and inputfile_B.off for booleans): a standard triangle OFF file.
• Edges (-e): an OFF-style file whose header is nv ne 0, followed by nv vertex lines x y z, then ne edge records 2 i j (the 2 is the number of endpoints, mirroring how OFF faces start with 3). See models/cube_subdiv_edges.off.
• Points (-p): an OFF-style file whose header is np 0 0, followed by np vertex lines x y z. See models/cube_subdiv_points.off.
• Existing tetrahedral mesh (for embedding a surface in it): a .tet or .msh (MEDIT) file passed as inputfile_B (ASCII). It is recognized by its extension, and a bool_opcode must still be given as a separator (it is ignored for embedding): mesh_generator surface.off U mesh.tet.

Output files

Inserting edges and points

Edges and points are forced into the tetrahedralization so they appear as genuine tet edges and vertices of the output. Internally, each input edge AB is pinned by two small triangles sharing AB, and each input point P by three small triangles meeting at P; these helper triangles are excluded from the surface tracking. Insertion implies -a (keep all cells), because the pinned edges/points do not bound a solid and would otherwise be discarded by the inside/outside classification.

# a surface plus extra edges and points, with exact tracking
mesh_generator models/sphere.off -t -r -a -e models/sphere_edges.off -p models/sphere_points.off

# only edges (no surface); the domain is the bounding box of the edges
mesh_generator -t -r -a -e models/cube_subdiv_edges.off

# only points
mesh_generator -t -r -a -p models/cube_subdiv_points.off

Tracking output back to the input

With -t the mesher writes sidecar files that map each input primitive to the output primitives it becomes. Use -a so the tracking is exact (every input primitive is realized). The three provenance sidecars share the same shape — keyed by the input primitive (coplanar triangle group / edge / point), each output primitive is given as an output tet index plus that tet's vertex indices (a face is tet v0 v1 v2, an edge is tet v0 v1, a point is tet vertex). All indices are 0-based into volume.tet's tet and vertex blocks; volume.tet.rational gives the exact coordinates of each output vertex in the same order.

volume.tet.groups — coplanar groups of the input triangles

The input triangles are partitioned into maximal edge-connected, exactly-coplanar groups (a flat region is one group; a triangle with no coplanar neighbour is a singleton group). Format:

# input_triangle_id coplanar_group_id
<num_input_triangles> <num_groups>
<input_triangle_id> <coplanar_group_id>
...

input_triangle_id is the triangle's index in the input OFF file (degenerate/collinear input triangles are dropped by the mesher, so this is not simply 0..n-1). A group is a 1-to-1 map to a single input triangle iff it has exactly one triangle.

volume.tet.triangleprov — triangle (face) provenance

One line per coplanar group, listing the output tet faces that tile it:

# group_id num_faces tet v0 v1 v2 ...
<num_groups>
<group_id> <num_faces> <tet v0 v1 v2> <tet v0 v1 v2> ...

Each face is given as an output tet index tet and the three output vertex indices v0 v1 v2 of that tet's face (v0 v1 v2 are three of tet tet's four vertices). A face on two exactly-coplanar surfaces that overlap on the same plane is listed under both groups. Summing the (distinct) face areas of a group reconstructs the group's area.

volume.tet.edgeprov — edge provenance

One line per input edge, listing the output tet edges that tile it (their union is exactly the input segment):

# edge_id num_out_edges tet v0 v1 ...
<num_input_edges>
<edge_id> <num_out_edges> <tet v0 v1> <tet v0 v1> ...

edge_id is the input edge's index in the -e file; each output edge is a tet index tet and the two output vertex indices v0 v1 of that tet's edge.

volume.tet.pointprov — point provenance

One line per input point:

# point_id tet out_vertex(-1 -1 if absent)
<num_input_points>
<point_id> <tet> <out_vertex>

out_vertex is the output vertex equal to the input point and tet an output tet containing it (both -1 if the point did not survive into the output — should not happen with -a).

Independent verifier

tests/verify_tracking.cpp re-derives all of the above in exact rational arithmetic and checks it (surface groups tile exactly, each edge is tiled by its output edges, each point equals its output vertex). It is a good reference for how to consume the sidecars:

mesh_generator models/cube_subdiv.off -t -r -a -e models/cube_subdiv_edges.off -p models/cube_subdiv_points.off
verify_tracking models/cube_subdiv.off volume.tet --strict \
    --edges models/cube_subdiv_edges.off --points models/cube_subdiv_points.off
# for an edges/points-only run (no surface), pass "none" as the surface:
verify_tracking none volume.tet --strict --edges models/cube_subdiv_edges.off

Examples on the bundled meshes

# 1) Volume mesh (polyhedral) of a closed surface
mesh_generator models/bunny.off                       # -> volume.msh

# 2) Tetrahedral mesh of the enclosed solid
mesh_generator models/bunny.off -t                    # -> volume.tet

# 3) Boolean union of two overlapping surfaces
mesh_generator models/cube_subdiv.off U models/sphere.off   # -> volume.msh

# 4) Exact surface tracking of a tetrahedralization
mesh_generator models/cube_subdiv.off -t -r -a        # -> volume.tet(.prov/.groups/.rational)

# 5) Insert edges and points, with tracking
mesh_generator models/sphere.off -t -r -a \
    -e models/sphere_edges.off -p models/sphere_points.off

# 6) Embed a surface into an existing tetrahedral mesh (.tet/.msh as inputfile_B).
#    The bool_opcode (here U) is a required separator and is ignored for embedding;
#    the surface must lie inside the tet mesh's domain.
mesh_generator models/Octocat.off U models/Octocat.bg.tet  # -> volume.msh (MEDIT: models/Octocat_medit.msh)

# 7) 2D: arrangement of a segment soup (see "2D mode" below)
mesh_generator models/2d/england.obj -2d              # -> triangulation.off
mesh_generator models/2d/siggraph_logo.obj -2d -r -v  # + exact rationals, verbose

2D mode

-2d runs a separate 2D pipeline: given a soup of segments in the plane — which may intersect, overlap, be duplicated or be degenerate — it computes their arrangement and outputs it as a triangulation in which every input segment is exactly a union of output triangle edges, plus the provenance mapping each input segment to the edges that tile it.

Internally: the input points are deduplicated and embedded in a Delaunay triangulation (a port of geogram's Delaunay2d onto this repository's exact predicate kernel), together with the four corners of the bounding box expanded by 10%. Each segment is then traced through the mesh; where it crosses an already-inserted segment, the intersection is created as an exact implicit point (segment-segment intersection) rather than a rounded coordinate. All geometric decisions go through the same exact predicates as the 3D path, so the output is likewise deterministic and byte-identical across platforms.

mesh_generator input.obj -2d [-r] [-v]

Input format

An OBJ line soup: v x y z vertex records and l i j line records (1-based indices; l may also carry a longer polyline, which is read as consecutive pairs). The z coordinate is ignored, so a planar curve exported from a 3D tool works directly. Two examples are bundled in models/2d/.

Duplicate points are merged, zero-length segments are dropped, and duplicate segments are processed once — all of them still appear in the provenance output.

Output files

Independent verifier

verify_tracking_2d re-derives the provenance from the geometry in exact rational arithmetic, sharing no code with the generator. For each input segment it checks that every output edge is exactly collinear with it and that the edges partition it with no gap and no overlap; it also checks that every input point maps to an output vertex with exactly equal coordinates, and that the triangulation is orientable and free of overlapping triangles.

mesh_generator models/2d/england.obj -2d -r
verify_tracking_2d models/2d/england.obj triangulation.off     # -> "TRACKING OK"

Viewing the result

tests/arrangement_to_vtu.py exports a run for ParaView as three files meant to be overlaid: the input segments, the output triangulation, and the output edges tagged with the id of the input segment they came from. Colouring the input and the output edges by the shared input_segment_id makes the provenance directly visible.

mesh_generator models/2d/siggraph_logo.obj -2d
python3 tests/arrangement_to_vtu.py models/2d/siggraph_logo.obj triangulation.off -o siggraph
# open siggraph_input.vtu, siggraph_output_tris.vtu, siggraph_output_edges.vtu in ParaView

Programmatic use

src/VolumeRemesher/2d/embed2d.h exposes embed_seg_in_tri_mesh, the 2D counterpart of embed_tri_in_poly_mesh: it takes flat arrays of segment coordinates and indices and returns the output triangles, exact rational vertex coordinates, and the segment and point provenance.

Programmatic API

From Python, see Python. From C++, the whole command line pipeline is driver.h (make_mesh, embed_mesh, arrange_2d, with the readers in io.h and the output snapshots and writers in mesh_output.h); main.cpp is a thin parser over it.

To embed a set of triangles (and, optionally, edges and points) into an existing tetrahedral mesh from C++, see embed.h (embed_tri_in_poly_mesh). It returns the polyhedral/tetrahedral output plus the same provenance information described above (faces on input, and per-edge / per-point provenance). Note: in the embedding path the edges/points — and the small helper-triangle apexes they generate — must lie inside the input tet mesh's domain (a debug assertion flags violations).

License

LGPL-3.0-or-later; see LICENSE. The default build, the Python wheels and the source distribution contain only LGPL-compatible components. The one exception is opt-in: third_party/graph_cut (the legacy -c graphcut classifier, see Inside/outside classification) is research-use-only and GPL code that is compiled only with -DVOLUMEREMESHER_WITH_GRAPHCUT=ON; a binary built that way must not be distributed under the LGPL.

Notes

• The pipeline uses exact indirect predicates, so the output is deterministic and byte-identical across platforms/compilers (this is checked in the integration tests), and identical between the command line tool and the Python package (checked by tests/python).

• We tested this code on macOS and Windows; it should also work on Linux.