Smooth Parametric Border Around a Nonplanar Mesh Boundary

[netemadi9]

Hey all!

I'd appreciate advice on the most suitable native MeshLib construction for this problem. I'll go over what is needed as well as the downstream uses of it to explain its constraints.

I start with an oriented triangle mesh and an ordered, closed 3D curve on its surface. After cutting out the selected region, I have a protected surface patch F whose boundary is C. I need to extend C outward into a configurable border B, ending at a new outer rim R. F and C must retain their geometry. The result F+B is initially an open surface with R as its boundary.

Border profile

Relative to a chosen construction axis d, the local profile has:

H — horizontal/transverse width: perpendicular to d, outward from the patch.
A — angled length: the actual inclined segment length, at angle θ above the transverse plane. Its transverse reach is A·cos(θ), and its axial rise is A·sin(θ).
V — vertical/axial height: parallel to d.

Optionally, an initial axial drop b before the horizontal segment, while retaining the original attachment at C. (This is not mandatory)

The problem

The boundary can be nonplanar and have concave turns, hairpins, narrow necks and nearby branches. It can also cross itself in projection along the chosen construction axis d while remaining separated in 3D (an overpass, not necessarily a true curve intersection.) Flattening it or welding spatially close branches loses information.

I'm not requesting an impossible constant-width, intersection-free offset for every curve. Local intersections between generated border triangles can be acceptable if their extent is small, the protected surface remains faithful, and downstream operations work. Broad folds, branch swaps, cracks and substantial unexplained width loss are not acceptable. When dimensions conflict with local clearance or curvature, I would like a principled, measurable compromise rather than a silent change of shape.

We currently use offsetContours, source mapping, constrained makeSpline smoothing, and guided MeshLib bridge/fill operations. Full H/A/V strips work well on ordinary boundaries; intermediate rows are retained where needed to preserve shape.

The remaining challenge is matching points between generated curves without distorting the border:

• Independent sampling creates mismatched vertex counts and triangle fans. Equal counts alone don’t guarantee correct correspondence.
• Forcing curves onto the outer rim’s samples distorted corners; combining every curve’s samples preserved outlines but created excessive tiny triangles.
• Other difficulties include extra contours, correspondence that doubles back, and coincident routed vertices. General welding can merge separate branches, while unconstrained smoothing can alter profile dimensions.

Is there a native workflow for coordinating these curves and their samples? The attachment C must remain fixed, and generated curves should retain the intended width, angle, height and smooth rim. Local transition fans and small, measured curve adjustments are acceptable.

The 3D contour offset operates in XY before restoring height. How should we handle overpasses... branches that cross in this flat view but remain separated in 3D... without accidentally joining them?

Lastly, we need this outer rim for two uses:

1. A minimum-wall-thickness exterior guide. This is a separate open surface on the intended outside of F, with a prescribed distance from F over the main region, ending exactly at R. It does not need to be fused to F+B. We now have a working approach for this guide. We need an unambiguous, ordered rim and a usable carrier surface. Tiny degenerate triangles, abrupt folds or branch ambiguity can still make the downstream construction harder.

2. Assembly with an outer shell. Eventually, an independently shaped outer shell must join F+B to form a closed object. For matching boundaries, exact seam stitching may be preferable; for overlapping surfaces or solids, a Boolean operation may be appropriate. We need consistent orientation, a manifold seam, no unintended gaps or degenerate triangles, and preservation of the protected geometry. I understand MeshLib's Boolean requirements concern self-intersections and inside/outside classification in the intersecting zone; accepting a local border crossing does not automatically mean it is admissible in that zone. A smooth rim alone is insufficient if its topology or correspondence is ambiguous.

A small example or recommended sequence of native operations would be extremely helpful. Thank you!

[netemadi9]

Oh, one more thing... how can we construct the generated curves together so that each horizontal/angled/vertical profile remains straight within its segment, while the border follows a winding boundary? We can allow small, controlled adjustments to the generated curves, but the attachment must remain fixed and the intended dimensions should remain faithful.

Here I illustrate our "goal" border as the first photo, with our attempt as the second photo. Note how our attempt introduces extra, non-uniform triangles and extra curvature to the horizontal, angled, and vertical profiles.

[CedricConday]

smoothRegionBoundary( Mesh&, const FaceBitSet& regionFaces, int numIters = 4 ) in MRMeshRelax.h does exactly this — it moves the region's boundary vertices to smooth the contour, with minor topology optimisation near the boundary. Use an even numIters; the solution oscillates back and forth.

[netemadi9]

Hey Cedric, thank you. It appears to have smoothed the existing boundary slightly, but it did not address constructing the border as described.

[CedricConday]

smoothRegionBoundary was the wrong primitive: it moves C, and you need C fixed. The native
construction for B is

MRMESH_API EdgeId extendHole( Mesh& mesh, EdgeId a, std::function<Vector3f(const Vector3f&)> getVertPos,
                              FaceBitSet* outNewFaces = nullptr );

in MRMeshFillHole.h. It adds a band of triangles from the current rim to new vertices placed at
getVertPos( existing vertex position ) and returns the edge of the new rim, so one call per
profile segment chains into the border: optional axial drop b, then the transverse run H, then
the inclined A (transverse A·cos(θ), axial A·sin(θ)), then the axial V. Each call moves the
rim by one straight step, which is what keeps every segment straight instead of accumulating the
curvature you see in the second photo, and C is never touched because only the newest rim is
extended. The result is left open at R, ready for whatever you cap it with.

For the direction: precompute one outward vector per boundary vertex from the loop tangent t and
your construction axis d, normalize( cross( t, d ) ) — take the loop from trackLeftBoundaryLoop
/ findLeftBoundary (MRRegionBoundary.h) and key the lookup by 3D position, since that is what
getVertPos receives. Nothing is projected along d anywhere in this, so an overpass or a hairpin
stays separated: two branches that collide in projection still have distinct positions and get
distinct directions.

[netemadi9]

Thank you! I will give this a try!

[Fedr]

The matching problem disappears if the curves are never created independently. Build the border row by row as mesh topology first, and assign positions afterwards: every row then has exactly the vertices of C, the correspondence is the mesh itself, and offsetContours - the only step in your pipeline that looks at the projection globally and can weld branches - is not used at all.

Recipe (all native, the attached script does it in ~100 lines):

1. findRightBoundary( mesh.topology ) gives the loop C with the hole to the left of every edge, vertices in order.
2. For every row of the profile call makeDegenerateBandAroundHole( mesh, e, &newFaces ). It returns the edge of the new hole opposite to e, so trackRightBoundaryLoop from that edge lists the new row in the same order as C. (extendHole( mesh, e, getVertPos ) does the same in one call when the new position depends only on the old one - that covers the drop b and the axial V segment, and extendHole( mesh, e, plane ) gives a planar rim - but H and A need the outward direction, so it is simpler to create degenerate rows and write the positions yourself.)
3. Per column i: transverse direction o_i = normal of the projected C (cross( d, tangent )), oriented into the hole - the sign follows from the loop orientation - and smoothed a few passes along the loop. Row positions are C_i + x*o_i + z*d with (x, z) the profile point of that row. Intermediate rows cost nothing.
4. Concave turns and narrow necks: measure instead of guessing. findSelfCollidingTriangles( MeshPart( mesh ) ), restricted to pairs touching the new faces, tells you exactly which columns collide. Scale their transverse reach by 0.85, re-place the rows (topology does not change), repeat; a Lipschitz clamp (reach may not change faster than 1:1 with arc length) keeps the rim smooth. This converges in ~10 iterations and gives you the achieved width per column as a number. The compromise is transverse only: axial heights, hence the rim height, are never changed.
   In a notch narrower than 2W (hairpin below) both sides stop at the middle and leave a slit. The true offset would fill the notch, but no 1:1 correspondence can express that; if you need it filled, weld the two rims where they meet (MeshBuilder::uniteCloseVertices( mesh, dist, /*uniteOnlyBd=*/true )) or skip the last shrink step and accept a small overlap.
5. Overpasses: nothing above looks at the projection globally and collisions are tested in 3D, so branches separated in 3D never interact. Helix example below: half a turn lies over the other half, 0 collisions, full width on every column.
6. Rim sampling: subdivideMesh with region = newFaces, subdivideBorder = true, a maxEdgeLen, and beforeEdgeSplit returning false for any edge with a face outside the band refines the stretched convex corners and R itself while F and C stay untouched (the F face set is identical afterwards).

Downstream

• The guide from Geometric Shell-Generation Problem #6885 needs only the ordered rim, and that is just the hole loop after the construction: one loop, |C| vertices, no ambiguity.
• Joining the outer shell: if the shell is open and its boundary runs near R, buildCylinderBetweenTwoHoles( mesh, eR, eShell, params ) (or stitchHoles) connects two loops with different vertex counts by a metric and gives a manifold seam without any resampling. If the shell overlaps F+B, make F+B closed first (extendHole/buildBottom plus fillHole on a far-away rim) and use boolean; check the band with findSelfCollidingTriangles before, and with the collision-driven reach there is nothing left to fix.

Results (profile A from your sketch plus b = 0.15; pip meshlib, a few seconds each)

Hairpin, legs 1.1 apart, reach H + A cos θ = 2.2, nonplanar C:

rows x columns	11 x 330, rim = 330 vertices, 1 loop
self-colliding pairs	919 -> 0 after 11 iterations
reach	full on 123 columns, reduced on 207, minimum 0.43 (= half the notch)

Helix strip, 1.5 turns, pitch 3.0 along d:

rows x columns	11 x 610, rim = 610 vertices, 1 loop
self-colliding pairs	0
reach	full on all 610 columns

Script: border6887.py (needs meshlib, numpy, matplotlib for the pictures).

[netemadi9]

Thank you so much!! I'm gonna give this a try - I love how it plays well with the geometric shell generation. Will keep you posted.

[netemadi9]

Fedr, thank you again - I think we're getting close to a solution.

Here's our pipeline at the moment, with a couple of remaining smart issues that I hope there is a simple solution for.

1. Cut and match boundaries. We use native isoline extraction and cutMesh to obtain the mesh opening on an outer shell, then cyclic arc-length correspondence to match it to the fixed rim. Our original preparation only used splitEdge, so it failed on a real case with 425 anatomy-boundary vertices versus 409 rim vertices.
2. Experimental unequal-count handling. We tested decimateMesh restricted to anatomy boundary edges, with endpoint-only collapses, no flips, and a 0.005 mm experimental simplification budget. We independently check ordered boundary-arc deviation, unchanged interior vertices/triangles, topology, and provenance, then split to the required rim count. Boundary adaptation passed its checks, with only 0.00373 mm maximum cut-boundary deviation. The remaining failure occurs during subsequent shaping and collision repair
3. Shape and weld. We use interpolateScalarsSmoothly with Unit weights to propagate boundary displacement, followed by shape-preserving Laplacian deformation with Cotan weights, graded positional attractors, an exactly pinned rim, and fixed upper anatomy outside the editable region. Finally, stitchContours creates an exact equal-count seam; no transition-band triangles are added.
4. Thickness-guided lift. A separate WeightedShell.meshShell surface provides thickness guidance. Signed nearest-surface queries and insertion-axis rays identify deficits. A custom small linear program combines nonnegative, compact geodesic displacement fields built using computeSurfaceDistances, with coherent cores, smooth falloff, rim fade, and a penalty for changing local anatomical relief. Movement is along the insertion axis. Candidates are screened for increased self/fitting intersections, creases, and triangle strain.
5. Independent fitting-clearance repair. We also detect actual anatomy–fitting/border crossings, calculate axial clearance demands, and translate coherent local cores with a harmonically blended transition using interpolateScalarsSmoothly. The rim stays fixed. Triangle degeneration/normal reversal and new self-intersections stop a candidate. Final validation checks closed/manifold topology, exact fitting/border preservation, and confirms collision candidates using exact stored-coordinate intersection tests.

Where it fails: Boundary adaptation and welding now succeed, but this same real case still fails geometrically. The thickness lift is rejected because it increases self-contact. The mesh we're lifting isn't perfectly smooth; it has peaks and valleys. The first fitting-repair candidate requires roughly 0.385 mm maximum axial clearance movement but reverses some triangle normals, so it stops before a second pass. Final validation finds 544 anatomy–fitting, 2 anatomy–border, and 1 self-intersection pair.
Three existing controls remain unchanged, and three deliberately oversampled-boundary controls pass. The experimental adaptation is not deployed.

What native MeshLib approach would you recommend for collision-aware, detail-preserving cervical deformation with an immutable rim...especially when axial lifting creates folds in the transition region? Would you favor local remeshing, a different constrained deformation formulation, or another boundary-adaptation strategy...?

[Fedr]

The failure you describe is structural, not a matter of tuning the lift, and I would stop deforming the fitting surface altogether. Two things are going on:

1. The axial direction is wrong exactly where you need it. On a wall tilted by α from the insertion axis, a normal deficit δ needs an axial move of δ / cos α; on an undercut it is unbounded, because the point has to travel past the bulge before the clearance appears. The cervical region is the steep region, so the "clearance demand" you compute there is huge, and lifting a whole wall by that amount is what folds the transition. On an undercut tooth with 2.4 mm of deficit the axial lift needs 9.3 mm and produces 1496 flipped triangles and 8576 self-intersections; the same vertices moved along the gradient of the distance field (towards / away from the nearest anatomy point) get there in 2.4 mm with 0 flips.
2. No displacement direction can fix the wedge. The fitting surface in the cervical wedge (wall + gum/plate) is a 2-D patch, and its target, the offset surface, has a crease there: the whole patch collapses onto a line. The gradient variant is fold-free but leaves 3232 pairs of overlapping slivers in that line; projecting onto a voxel guide leaves 1296. Detail-preserving Laplacian shaping, attractors and LP-weighted fields all move vertices; the wedge needs different triangles, so any deformation formulation ends here.

Formulate the constraint as a set operation instead. "At least t from the anatomy" means

fitting surface  =  upper boundary of   ( solid under D )  ∪  ( points within t of the anatomy )

and that is a plain boolean( Union ) on two closed meshes:

• Solid under D: extendHole( D, rimEdge, Plane3f( d, far ) ) then fillHole on the new hole - a vertical skirt from R down to a far plane plus a cap. R is untouched by construction. Collect the skirt and cap faces with outNewFaces.
• Clearance solid: offsetMesh( MeshPart( anatomy ), t, OffsetParameters{ voxelSize } ). It is closed by construction; for an open anatomy it is the two-sided shell, and the lower sheet lies inside the D-solid and disappears in the union. Build it from the anatomy alone (not from F+border), so it stays farther than the border width from R; two solids sharing the rim would give coincident skirts, which a boolean cannot handle. Your WeightedShell guide works too if it is closed the same way and does not touch R.
• boolean( solidD, clearanceSolid, BooleanOperation::Union, nullptr, &mapper ), then deleteFaces( mapper.map( skirtAndCap, MapObject::A ) ) and pack().

Everything you were validating by hand comes out of the construction: the boundary is R to 1e-8, clearance ≥ t up to the voxel error, no self-intersections, the anatomy-fitting crossings are gone (a crossing is just a negative deficit), and D is bit-identical wherever it already had clearance - 65-74 % of its vertices in the tests, the rest replaced by the offset surface with the boolean's own seam. No transition band, no LP, no candidate screening; both tests below run in 7 s together. If the crease where D meets the offset must be soft, subdivide around the intersection contour and positionVertsSmoothly only there, with the offset taken at t + margin so the smoothing never dips under t.

Tests (pip meshlib, scripts attached)

Undercut tooth of revolution on a plate, fitting shell at 1 mm where 3 mm is required, rim = plate boundary fixed:

repair	anatomy-fitting crossings	self-intersections	flipped triangles	min clearance	max move
axial lift + harmonic blend, 6 rings	0	8576	1496	3.00	9.3
gradient displacement onto the offset + smooth ring	0	3232	0	2.33	2.4
nearest point on a voxel guide + smooth ring	0	1296	254	2.97	2.6
union of solids	0	0	-	2.99	rim deviation 4e-15, 74 % of D untouched

Dipped design shell over the #6885 anatomy (moderate slope, 38°, dip 4.5 mm reaching 1.5 mm into F): here the axial lift does work (0 crossings, 7 flips), the union gives 0 crossings, 65 % of D untouched, and the 11 self-intersecting pairs it reports are the seam slivers of the input guide itself.

On your steps 1-2: with the union the fitting surface's boundary is R, so nothing has to be resampled or decimated to a matching count before stitchContours; and if you ever need to connect two loops with different counts, buildCylinderBetweenTwoHoles / stitchHoles do it by a metric without touching either loop.

Scripts: tooth6887.py (steep test, three deformation variants), lift6887.py (the #6885 data with a dip), union6887.py (the set-operation version on both). clearance6887_prototype.zip

[netemadi9]

Fedr, thank you. The union makes sense, and your point about lifting along the axis on steep walls matches what we saw. Since first writing this, we ran your prototype and several alternatives on some cases, so I've narrowed my questions.

What we learned

• The crossings in our failing case were not at the rim. They were the occlusal surface of a tooth placed too low, 1–5.7 mm from R. Our seating made it worse: it bent a 3.5–6 mm band of the tooth onto R, which pulled the occlusal surface down about 0.6 mm. We are shortening it so it bends only the neck.
• The "more edge points than rim points" failure came from our cut, not from the pairing. We measured each vertex's height above the nearest rim segment, and on a nearly level cervical shelf the 0.6 mm isoline wandered: about a third longer than R (583 points vs 513), plus a 6 µm island. A smooth rim-height field gives one loop with well under |R| points.
• Your union reproduced your numbers on our MeshLib (3.1.4) and removed every crossing on the real case. On real data it needed three additions: restoring the kept vertices' coordinates through old2newVerts (the boolean moved some by about 1e-8 mm, which broke our exact rim weld), removing zero-area faces away from the boundary, and getPlaneFillMetric for the skirt cap (the default metric folded it). It does fill fossae with the offset surface, so we keep it as a last-resort safety net.

1. The lift
Its purpose is a targeted increase of the occlusal anatomy where the tooth is thinner than the minimum-thickness surface E. Our rules:

• it only raises, never lowers (a tooth placed too high stays high);
• only the thin area moves, plus a small blend: if 2 mm of the occlusal table is thin, only that area lifts;
• inside that area the relief is kept: fossae stay deep and may flatten only slightly;
• crossings with the fitting surface must be cleared; E is the target.

Our prototype is a thin-plate displacement s along the insertion axis with s ≥ (deficit to E + 0.03 mm), solved by an active set in scipy. It clears everything, but it lifted nearly the whole table (median 0.6 mm), and the deepest fossae lost 24–30% of their depth. We are now restricting it to the thin region plus a blend.

Is there a native MeshLib way to solve this: a smooth displacement over a vertex region, with lower bounds (s ≥ d) on part of it and zero value and slope at its edge? And how would you keep the fossae from filling?

2. The tooth's edge and the rim
With the smooth cut, the cut edge now has fewer points than the rim, so our one-to-one pairing works. Two questions:

• Is there a MeshLib way to extend R's heights smoothly over the plane across the axis (for example a harmonic extension), so the cut line stays a fixed height above the local rim without wandering?
• We also tried stitchHolesNicely for joining loops with different counts. The join works, but its subdivision split the R and K edges. When we blocked those splits with beforeEdgeSplit, subdivideMesh never returned in our runs if a blocked edge was longer than maxEdgeLen. Setting maxEdgeLen to at least the longest loop edge avoided it. Is that expected?

[Fedr]

Good to hear the union holds on real data; getPlaneFillMetric for the cap is the right call (a planar hole with a non-convex outline can fold under the default metric).

1. The lift

The native solver for "smooth displacement over a region, handles fixed, zero value and zero slope at the edge, relief kept" is Laplacian (MRLaplacian.h):

Laplacian lap( mesh );
lap.init( freeVerts, EdgeWeights::Cotan, VertexMass::Unit, RememberShape::Yes );
for ( auto v : handles ) lap.fixVertex( v, mesh.points[v] + s[v] * d );
lap.apply();

With RememberShape::Yes the free vertices are solved so that their Laplacian coordinates stay what they were, which is exactly detail preservation; the system is separable per coordinate, so handles moved along d move everything along d only (tangential drift 0 in the test). The lower bound s ≥ deficit is an active set on the same object: after apply(), any free vertex below E + margin becomes a handle at E + margin (fixVertex again, apply again); it converged in one or two rounds here. applyToScalar does the same for a scalar field if you prefer to solve for s itself. One caveat: a bi-harmonic blend overshoots like a spline (the tests show 0.05-0.12 mm dips next to the lifted zone), so "never lowers" is not automatic; if you need the maximum principle, interpolateScalarsSmoothly (harmonic) never leaves the range of the fixed values, at the price of a C0 kink at the region edge.

The fossae are not a solver problem. Depth is z(cusp) - z(floor), the floor must reach E + margin, so either the cusps rise by the deficit or the depth shrinks by it; your 24-30 % is that deficit divided by that depth. The only free choice is where the lift is spread: give every vertex within a radius ρ of the thin area the maximum deficit found in its ρ-neighbourhood (an envelope), and let the Laplacian blend the edge. ρ below the fossa size fills the fossa; ρ reaching the cusps keeps the depth exactly and raises the whole table around it. Nothing in between keeps both.

Occlusal table with four cusps and a central fossa 0.5 mm below E, margin 0.03, 6 blend rings:

handles	table moved	fossa depth	min (z - E) after	dip below original
per-vertex deficit (thin-plate over the thin area)	7 %	2.38 -> 1.85 (-22 %)	+0.03	-0.06
envelope ρ = 2 mm	31 %	2.38 -> 1.83 (-23 %)	+0.03	-0.05
envelope ρ = 4.5 mm (reaches the cusps)	84 %	2.38 -> 2.38 (0 %)	+0.03	-0.12

2. Cut line at a fixed height above the rim

Harmonic extension is one call: fix the rim vertices at their height in a VertScalars, set InterpolateScalarsParams::region to all other vertices with EdgeWeights::Cotan, and interpolateScalarsSmoothly. The cut is then extractIsolines( topology, z - h, 0.6 ). On the #6885 guide (rim 1344 vertices) this gives one loop with 856 points at 0.6 (901 at 0.3, 817 at 1.0), no islands. If "height above the rim" was a stand-in for "distance from the rim", computeSurfaceDistances from the rim vertices and an iso-distance cut is the other one-liner.

3. subdivideMesh not returning with protected edges

Expected. Subdivision stops when every edge in the region is at most maxEdgeLen. A triangle on a protected edge of length L has another side of at least L/2, so with maxEdgeLen < L that state is unreachable, only maxEdgeSplits ends the loop, and every split re-queues the rejected edge, so the time grows superlinearly: in a small repro (loops of 60 and 40 vertices, protected loop edges 0.5-0.6 mm, maxEdgeLen 0.3) 3000 splits take 0.5 s, 20 000 splits do not finish in a minute, 100 000 run out of memory. Your rule is the right one: maxEdgeLen not below the longest protected edge. If you need a finer patch than that, split R and K yourself first, identically on both sides of the weld, and protect the result.

Scripts: lift_occlusal.py, rimfield.py, stitch_hang.py. lift6887_prototype.zip