Skip to content

Isotropic remeshing

isotropic_remeshing remeshes a triangle mesh towards a uniform target edge length (Botsch and Kobbelt 2004). Each iteration splits the edges longer than 4/3 of the target, collapses those shorter than 4/5 of it, swaps edges to bring the valences closer to 6 (4 on the boundary), and moves every vertex towards the average of its neighbours, within its tangent plane. Nothing holds the surface by default: with an envelope, no operation may move it farther than a given distance from the input; with features, curves and vertices the remeshing keeps.

Use it for uniform, well-shaped triangles at a chosen edge length, finer or coarser than the input. To reduce a mesh to a vertex budget while keeping its shape, use QSlim; for a volume mesh, TetWild.

A first run

The bunny (35,947 vertices in the file, 34,834 in its triangles) remeshed to edges of 1% of its bounding box diagonal:

{
  "application": "isotropic_remeshing",
  "input": "bunny.obj",
  "output": "bunny_remeshed.obj",
  "length_rel": 0.01,
  "report": "bunny_report.json"
}

wmtk_app -j isotropic_remeshing.json

The input is read relative to the file; the remeshed mesh, bunny_remeshed.obj, and the report, bunny_report.json, are written to the working directory.

bunny = wm.read_mesh(str(DATA / "bunny.obj"))
outputs, report = wm.isotropic_remeshing({"input": "bunny", "length_rel": 0.01}, {"bunny": bunny})
remeshed = outputs["mesh"]
print(len(bunny["V"]), "->", len(remeshed["V"]), "vertices")
print(f"target length {report['target_length']:.5f}, average {report['avg_length']:.5f}")

On files instead, with the same options as the JSON file:

wm.isotropic_remeshing({
    "input": str(DATA / "bunny.obj"),
    "output": "bunny_remeshed.obj",
    "length_rel": 0.01,
    "report": "bunny_report.json",
})

wmtk::MeshMap meshes{{"bunny", wmtk::io::read_mesh(data / "bunny.obj")}};
auto [outputs, report] =
    isotropic_remeshing({{"input", "bunny"}, {"length_rel", 0.01}}, meshes);
const wmtk::MeshData& remeshed = outputs.at("mesh");
std::cout << meshes.at("bunny").V.rows() << " -> " << remeshed.V.rows() << " vertices\n";
wmtk::io::write_mesh("bunny_remeshed.msh", remeshed);

Link wmtk::isotropic_remeshing.

It prints 35947 -> 7929 vertices. The run is sequential, so the result is the same every time.

The input on the left, the result on the right: drag to rotate, scroll to zoom. Every example is in the gallery.

Key settings

The target length: length_rel, length_abs

length_abs is the target edge length; when it is negative (the default), length_rel gives it relative to the bounding box diagonal; when both are negative, the target is the input's average edge length. Since edges between 4/5 and 4/3 of the target are left alone, the average edge of the result is somewhat longer than the target. On the bunny:

for length_rel in [0.005, 0.01, 0.02]:
    outputs, report = wm.isotropic_remeshing(
        {"input": "bunny", "length_rel": length_rel}, {"bunny": bunny})
    print(f"length_rel {length_rel:g}: {report['#v']} vertices, average edge "
          f"{report['avg_length'] / report['target_length']:.2f} x the target")
length_rel vertices average edge / target
0.005 30,313 1.18
0.01 7,929 1.15
0.02 2,034 1.12

Halving the length multiplies the vertex count by about four, as it should on a surface.

The iterations: max_iterations

max_iterations (default 3) is the number of iterations of splits, collapses, swaps and smoothing; all of them run. The first does most of the work. On the bunny, the share of triangles with an angle under 30 degrees:

print("input: smallest angle %.1f, %.2f%% of the triangles under 30 degrees" % quality(bunny))
for max_iterations in [1, 3, 10]:
    outputs, _ = wm.isotropic_remeshing(
        {"input": "bunny", "length_rel": 0.01, "max_iterations": max_iterations},
        {"bunny": bunny})
    print(f"max_iterations {max_iterations}: smallest angle %.2f, %.2f%% under 30 degrees"
          % quality(outputs["mesh"]))
under 30 degrees
input 6.77%
1 iteration 1.66%
3 0.76%
10 0.69%

The smallest angle stays at 0.02 degrees: the worst triangles are along the frozen boundary, and more iterations do not reach them (see the boundary).

The envelope: eps_rel

eps_rel is the envelope's thickness, relative to the input's bounding box diagonal. The default, -1, is no envelope: smoothing moves every vertex within its tangent plane, which on a curved surface or at a sharp edge is off the surface, and the iterations add up. A box (piece_0.obj, 8 vertices), remeshed with length_rel 0.05, keeps none of its corners: they round off. With an envelope, an operation that would leave part of the surface outside it is refused:

for eps_rel in [-1, 0.01, 0.001]:
    outputs, _ = wm.isotropic_remeshing(
        {"input": "box", "length_rel": 0.05, "eps_rel": eps_rel}, {"box": box})
    print(f"eps_rel {eps_rel:g}: farthest vertex {farthest_from_box(outputs['mesh']):.2g} "
          "from the box")
eps_rel farthest vertex from the box
-1 (none) 0.023
0.01 0.0066
0.001 0.00069

The diagonal is 1.22, so each envelope holds with room to spare. To keep sharp edges sharp, rather than within a distance, make them features.

The envelope test is exact by default, and costs most of the run: on the bunny, an envelope of 0.2% makes the first run close to a hundred times slower, nearly all of it in the smoothing, and the mesh stays somewhat finer where collapses are refused. use_sample_envelope tests against the same envelope by sampling: there, about thirty times faster than the exact test, but it refuses more operations and leaves more vertices. With more than one thread (num_threads), the result changes slightly from run to run (7,878 and 7,916 vertices for the first run with 4 threads, 7,929 sequentially); 0, the default, runs sequentially.

The open boundary: freeze_boundary

By default (freeze_boundary true), the vertices of an open boundary are frozen: they never move, and no collapse touches them. Splits still refine the boundary edges. The bunny has 223 boundary edges, around five holes:

print("input:", boundary_edges(bunny), "boundary edges")
for options in [{}, {"freeze_boundary": False}, {"boundary_eps_rel": 0.001}]:
    outputs, _ = wm.isotropic_remeshing(
        {"input": "bunny", "length_rel": 0.01, **options}, {"bunny": bunny})
    print(options, boundary_edges(outputs["mesh"]), "boundary edges, smallest angle %.2f" %
          quality(outputs["mesh"])[0])
boundary edges smallest angle
frozen (default) 230 0.02
freeze_boundary false 102 1.94
boundary_eps_rel 0.001 196 0.00

Frozen, the boundary keeps every input vertex, and the triangles next to it keep their short edges: they are the thinnest of the result. With freeze_boundary false, the boundary is remeshed like the rest, but nothing holds it: the holes can change shape. A positive boundary_eps_rel remeshes the boundary inside a tube around the input's boundary edges. In all three runs, the thinnest triangles are along the boundary: the swaps by valence leave them there, and swap_quality removes most of them.

Features: features, feature_fields

Features are curves and vertices the remeshing keeps. A feature edge is never swapped; a split leaves two feature edges; a vertex on a curve moves only along it, and collapses only along it. The ends of a curve, the vertices where three or more curves meet or where the edge id changes, and fixed vertices stay where they are. They come in two forms:

  • features lists them: on the command line, a JSON file {"edges": [[a, b], ...], "edge_ids": [id, ...], "vertices": [[v, id], ...]}, its indices vertices of the input file, each edge an edge of the input; in memory, a mesh on the input's vertex rows, whose edges are the curves (with an Int64 cell field edge_id) and whose points are the fixed vertices (with vertex_id). The ids are carried to the output: both halves of a split edge keep theirs. With an envelope, a tube of the same thickness holds the curves.
  • feature_fields names per-face fields, labels, whose boundaries become curves: on the command line and on files, fields of a .msh input; in memory, triangle attributes of the input. Every output triangle has the value of the input triangle it comes from (an Int64 label has to be within 2^53 in magnitude, which the remeshing carries exactly), and a degenerate input triangle, which the reader drops, takes its value with it. With an envelope, a tube of the same thickness holds these curves too.

The box again, without an envelope, its 12 sharp edges found by the angle between their two triangles, each with an id, and its 8 corners fixed:

def sharp_edges(mesh, degrees=30):
    """The edges of a closed manifold mesh whose two triangles' normals are more than `degrees`
    apart."""
    V, F = mesh["V"], mesh["cells"]["triangle"]["connectivity"]
    N = np.cross(V[F[:, 1]] - V[F[:, 0]], V[F[:, 2]] - V[F[:, 0]])
    N /= np.linalg.norm(N, axis=1, keepdims=True)
    E = np.sort(F[:, [0, 1, 1, 2, 2, 0]].reshape(-1, 2), axis=1)
    face = np.repeat(np.arange(len(F)), 3)
    order = np.lexsort((E[:, 1], E[:, 0]))
    E, face = E[order], face[order]
    # Each edge is there twice, in consecutive rows: once for each of its triangles.
    edges, f0, f1 = E[0::2], face[0::2], face[1::2]
    return edges[(N[f0] * N[f1]).sum(1) < np.cos(np.radians(degrees))]


edges = sharp_edges(box)
# The corners: the vertices of three or more sharp edges.
corners = np.flatnonzero(np.bincount(edges.ravel(), minlength=len(box["V"])) >= 3)

In memory, the curves are a mesh:

# The curves as a mesh on the vertex rows of the box: its edges, each with an id, and its
# points, the fixed vertices, each with an id too.
curves = {
    "V": box["V"],
    "cells": {
        "edge": {"connectivity": edges.astype(np.int64),
                 "attributes": {"edge_id": {"values": np.arange(len(edges), dtype=np.int64)}}},
        "point": {"connectivity": corners.astype(np.int64).reshape(-1, 1),
                  "attributes": {"vertex_id": {"values": corners.astype(np.int64)}}},
    },
}
outputs, _ = wm.isotropic_remeshing(
    {"input": "box", "length_rel": 0.05, "features": "curves"}, {"box": box, "curves": curves})
mesh, feature_edges = outputs["mesh"], outputs["edges"]
edge_id = feature_edges["cells"]["edge"]["attributes"]["edge_id"]["values"]
vertex_id = mesh["vertex_attributes"]["vertex_id"]["values"]
print(f"{len(edges)} sharp edges -> {len(edge_id)} feature edges, "
      f"{len(np.unique(edge_id))} ids, {(vertex_id >= 0).sum()} fixed vertices")
print(f"farthest vertex {farthest_from_box(mesh):.2g} from the box")

On files, a features file, and a .vtu output for the curves and the ids:

# The same curves as a features file: indices are vertices of the input file.
with open("box_features.json", "w") as f:
    json.dump({"edges": edges.tolist(),
               "edge_ids": list(range(len(edges))),
               "vertices": [[int(v), int(v)] for v in corners]}, f)
wm.isotropic_remeshing({
    "input": str(DATA / "piece_0.obj"),
    "output": "box_remeshed.vtu",  # and box_remeshed_edge.vtu, the feature edges
    "length_rel": 0.05,
    "features": "box_features.json",
})

wmtk::MeshData box = wmtk::io::read_mesh(data / "piece_0.obj");
// The sharp edges: those whose two triangles' normals are more than 30 degrees apart.
const auto& F = box.cell(wmtk::CellType::Triangle).connectivity;
std::map<std::pair<int64_t, int64_t>, std::vector<Eigen::Vector3d>> normals;
for (Eigen::Index f = 0; f < F.rows(); ++f) {
    const Eigen::Vector3d a = box.V.row(F(f, 0)), b = box.V.row(F(f, 1)),
                          c = box.V.row(F(f, 2));
    for (int k = 0; k < 3; ++k) {
        const int64_t u = F(f, k), v = F(f, (k + 1) % 3);
        normals[{std::min(u, v), std::max(u, v)}].push_back((b - a).cross(c - a).normalized());
    }
}
std::vector<std::pair<int64_t, int64_t>> sharp;
std::vector<int> degree(box.V.rows(), 0);
for (const auto& [e, n] : normals) {
    if (n.size() == 2 && n[0].dot(n[1]) < std::cos(30.0 / 180.0 * 3.14159265358979)) {
        sharp.push_back(e);
        ++degree[e.first];
        ++degree[e.second];
    }
}
// The curves as a mesh on the vertex rows of the box: its edges, each with an id, and its
// points, the corners (on three or more sharp edges), each with an id too.
wmtk::MeshData curves;
curves.V = box.V;
wmtk::CellData& edges = curves.cells[wmtk::CellType::Edge];
edges.connectivity.resize(sharp.size(), 2);
wmtk::DataMatrix<int64_t> edge_id(sharp.size(), 1);
for (size_t i = 0; i < sharp.size(); ++i) {
    edges.connectivity.row(i) << sharp[i].first, sharp[i].second;
    edge_id(i, 0) = int64_t(i);
}
edges.attributes["edge_id"].values = edge_id;
std::vector<int64_t> corners;
for (size_t v = 0; v < degree.size(); ++v) {
    if (degree[v] >= 3) corners.push_back(int64_t(v));
}
wmtk::CellData& points = curves.cells[wmtk::CellType::Point];
points.connectivity.resize(corners.size(), 1);
for (size_t i = 0; i < corners.size(); ++i) points.connectivity(i, 0) = corners[i];
points.attributes["vertex_id"].values = points.connectivity; // a corner's id: its index

const auto [box_outputs, box_report] = isotropic_remeshing(
    {{"input", "box"}, {"length_rel", 0.05}, {"features", "curves"}},
    {{"box", box}, {"curves", curves}});
std::cout << sharp.size() << " sharp edges -> "
          << box_outputs.at("edges").cell(wmtk::CellType::Edge).connectivity.rows()
          << " feature edges\n";

It prints 12 sharp edges -> 128 feature edges, 12 ids, 8 fixed vertices, and every vertex is on the box: the farthest is at distance 0. On the command line, "features": "box_features.json" names the same file.

The same 12 edges are the boundaries of a label per side of the box: written in a .msh file as a per-face field, feature_fields finds them, and the output has 128 feature edges again, each triangle with its side's label:

# The same curves as the boundaries of a label per face of the box, in a .msh file.
V, F = box["V"], box["cells"]["triangle"]["connectivity"]
normal = np.cross(V[F[:, 1]] - V[F[:, 0]], V[F[:, 2]] - V[F[:, 0]])
axis = np.abs(normal).argmax(1)
patch_id = 2 * axis + (normal[np.arange(len(F)), axis] > 0)  # 0 to 5: the side it faces
wm.write_mesh("box_labels.msh", {"V": V, "cells": {"triangle": {
    "connectivity": F, "attributes": {"patch_id": {"values": patch_id.astype(float)}}}}})
wm.isotropic_remeshing({
    "input": "box_labels.msh",
    "output": "box_labels_remeshed.vtu",  # with patch_id on every triangle
    "length_rel": 0.05,
    "feature_fields": ["patch_id"],
})

Both write a .vtu: the triangles, with the per-face fields and the point field vertex_id (the fixed vertices' ids, -1 elsewhere), and <stem>_edge.vtu next to it, the feature edges with their edge_id (-1 for a boundary of feature_fields). Any other format holds the triangles only, with a warning.

The swaps: swap_quality

By default an edge is swapped where that brings the valences of its four vertices closer to 6. That leaves alone a thin triangle whose swap would not change the valences much: next to a frozen boundary or a feature curve, where the curve itself cannot be swapped. With swap_quality, an edge is swapped where that raises the worse shape quality of its two triangles. On the bunny:

for options in [{}, {"swap_quality": True}, {"swap_quality": True, "boundary_eps_rel": 0.001}]:
    outputs, _ = wm.isotropic_remeshing(
        {"input": "bunny", "length_rel": 0.01, **options}, {"bunny": bunny})
    print(options, "smallest angle %.2f, %.2f%% under 30 degrees" % quality(outputs["mesh"]))
smallest angle under 30 degrees
valence (default) 0.02 0.76%
swap_quality 7.41 0.49%
swap_quality, boundary_eps_rel 0.001 14.84 0.42%

The valences are then less regular. Turn it on when the smallest angle matters more than the valences, and with features.

use_link_condition (default true) refuses a collapse that would change the topology: the genus, the boundary loops, or the manifoldness. Turned off, a coarse target can change the surface itself: sphere.obj remeshed with edges of half its diagonal ends as a single triangle instead of a tetrahedron. Turn it off only for inputs whose topology is wrong already: on double_sphere.obj, two overlapping spheres joined by non-manifold edges, the link condition refuses many collapses at a coarse target, and leaves more vertices and thinner triangles than without it. Swaps are refused on non-manifold edges either way.

Outputs

  • The remeshed mesh: output on the command line (default out.obj), outputs["mesh"] in memory. With features, a .vtu output and its _edge.vtu, and in memory outputs["mesh"] with the vertex field vertex_id and the per-face fields of feature_fields (each in the type it was given in, with its labels), and outputs["edges"], the feature edges on the same vertex rows with edge_id. A result goes back in as it is: mesh as the input, edges as the features. One made with feature_fields goes back in with its per-face fields named in feature_fields again (mesh holds no uv or normal): an input may carry no triangle field that feature_fields does not name.
  • With report, a JSON file with the target length, the average, smallest and largest edge length, the valences, the vertex and triangle counts, the non-manifold edges and vertices, the number of operations refused because they would fold the surface, and the run time. In memory, the report is returned.
  • With history_output, the history of the operations, from which tracking maps points between the input and the output.

Recipes and pitfalls

A label map

A mesh with a label per triangle, remeshed with the boundaries between labels kept as curves: topo_annots.msh, a flat square of 2,601 vertices with six labels in the field tag_0. The square's outline is in a tube rather than frozen, so that it is remeshed to the target too, and the swaps go by quality, since the label boundaries are curves:

{
  "application": "isotropic_remeshing",
  "input": "topo_annots.msh",
  "output": "topo_annots_remeshed.vtu",
  "feature_fields": ["tag_0"],
  "length_abs": 2,
  "boundary_eps_rel": 0.001,
  "swap_quality": true
}

wmtk_app -j isotropic_remeshing_labels.json

wm.isotropic_remeshing({
    "input": str(DATA / "topo_annots.msh"),
    "output": "topo_annots_remeshed.vtu",
    "feature_fields": ["tag_0"],
    "length_abs": 2,
    "boundary_eps_rel": 0.001,
    "swap_quality": True,
})

The 310 edges between labels become curves; the result has 1,177 triangles, each with its label in topo_annots_remeshed.vtu, and 92 feature edges in topo_annots_remeshed_edge.vtu.

Pitfalls

  • A misspelled option is ignored. isotropic_remeshing runs with the default of an option it does not know; wmtk_app and the Python bindings log a warning (isotropic_remeshing has no option ...). Editors catch it beforehand with the JSON Schemas.
  • The surface drifts without an envelope. Smoothing is only tangential to the current mesh: set eps_rel when the result has to stay close to the input, and features for edges that must stay sharp.
  • The target length is not reached with an envelope. A collapse that would leave the envelope is refused, so where the envelope is thin compared with the target length, the mesh stays finer, and its triangles worse. Thicken the envelope or accept the finer mesh.
  • Thin triangles next to curves or a frozen boundary. Turn on swap_quality; and cleanup_angle, which repairs the triangles under a given angle at the end of every iteration.
  • The features are lost. Only a .vtu output holds the feature edges, their ids and the per-face fields; an .obj holds the triangles, and the run warns.
  • A labelled result refused as the next input. In memory, outputs["mesh"] carries the per-face fields of feature_fields, and an input may carry no triangle field that feature_fields does not name (Unknown mesh attribute): pass the same per-face fields again. The result holds no uv or normal, only their seams, in edges.
  • Inconsistently oriented input. The fold check (check_inversions) counts every flip of orientation inside a region of the feature curves as a fold -- not one along a curve between two regions -- so most operations there are refused. Orient the input first.
  • The output is not next to the JSON file. wmtk_app reads input (and features) relative to the JSON file, but writes output and report relative to the working directory.

See also