Skip to content

smoothing

triwarp.smoothing

Laplacian smoothing filters, for vertex positions and for per-vertex scalar fields.

Most of the module moves geometry: filter_laplacian, filter_taubin, filter_humphrey, filter_neighborhood_average, filter_mut_dif_laplacian and filter_implicit_fairing all diffuse vertex positions through the same row-stochastic 1-ring operator, differing in the time integration and in what they do to counteract shrinkage. smooth_region and its sharp-boundary variant instead solve a Dirichlet problem over a region, holding the rest of the mesh fixed.

Three functions break that pattern by working on the normal field instead of positions, which is what lets them keep a crease sharp: filter_normals diffuses face normals with a crease gate, filter_two_step then refits the vertices to them, and filter_sharpen runs the whole idea backwards to sharpen.

A second group relaxes toward something other than a Laplacian residual, which is what lets each member fix a failure the filters above cannot see: equalize_triangle_areas evens out triangle areas, relax_keep_volume removes the shrinkage locally instead of rescaling it away at the end, relax_approx fits a plane or a quadric to a whole geodesic neighbourhood rather than averaging a 1-ring, and filter_spikes moves only the vertices that fail an angle test. The first three take a vertex region and a max_displacement bound, so a relaxation can be confined to where it is wanted and kept within a tolerance of the surface it started from; filter_spikes needs neither, because the set it touches is the answer to its own test. smooth_region_boundary completes the region trio by smoothing the region's rim curve, where the two above it smooth across the rim or inside it.

The verb tracks the mechanism, not the group, which is why that second group spans three of them. filter_* runs a fixed operator to a schedule -- an assembled Laplacian, a normal-field pass, a windowed Taubin pair -- so the answer is a function of the operator and the iteration count. relax_* iterates against a geometric objective and re-derives its target every pass, which is why those take a max_displacement bound: there is nothing in the mechanism that keeps the result near the input. equalize_triangle_areas and smooth_region* name their objective outright because there is only one of each. By that axis filter_spikes is correctly a filter_*: it runs the same fixed 1-ring operator the filters above it do, and only its selection is geometric. It sits in the paragraph above because its selection is what makes it useful, not because it is a relaxation.

filter_scalar_laplacian runs the same operator over a per-vertex scalar field rather than positions. Capping how fast such a field may vary along an edge — the other half of turning a raw scalar into a usable sizing field — is not a smoothing filter at all but a one-sided Lipschitz projection, and lives in shortest_path_envelope.

filter_laplacian 

filter_laplacian(
    vertices: array[vec3],
    faces: array[int32],
    lamb: float = 0.5,
    iterations: int = 10,
    implicit_time_integration: bool = False,
    volume_constraint: bool = True,
    laplacian_operator: BsrMatrix[float32] | None = None,
) -> wp.array[wp.vec3]

Laplacian mesh smoothing (Vollmer et al.; Desbrun et al. implicit fairing).

Diffuses each vertex toward the mean of its 1-ring neighbors. With implicit_time_integration the diffusion is solved implicitly (backward Euler, unconditionally stable); otherwise an explicit forward-Euler step is used (stable for lamb <= 1). Mirrors trimesh.smoothing.filter_laplacian but returns a new vertex array instead of mutating the mesh in place.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
lamb float

Diffusion speed. 0 leaves the mesh unchanged; larger values smooth faster.

 0.5
iterations int

Number of smoothing passes.

 10
implicit_time_integration bool

If True solve ((1 + lamb) I - lamb L) V' = V each pass via conjugate gradient. If False apply the explicit step V' = V + lamb (L V - V).

 False
volume_constraint bool

If True rescale the mesh after each pass to preserve its initial volume, counteracting Laplacian shrinkage.

 True
laplacian_operator BsrMatrix[float32] | None

Optional precomputed row-stochastic operator (see laplacian). Autogenerated (uniform weights) when None.

 None

Returns:

Type Description
array[vec3]

Smoothed (n_vertices,) vertex positions on vertices.device.

Raises:

Type Description
RuntimeError

If vertices and faces are not all on one device.
See Also

filter_taubin filter_humphrey filter_implicit_fairing trimesh.smoothing.filter_laplacian

Source code in triwarp/smoothing.py 
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def filter_laplacian(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    lamb: float = 0.5,
    iterations: int = 10,
    implicit_time_integration: bool = False,
    volume_constraint: bool = True,
    laplacian_operator: wps.BsrMatrix[wp.float32] | None = None,
) -> wp.array[wp.vec3]:
    """
    Laplacian mesh smoothing (Vollmer et al.; Desbrun et al. implicit fairing).

    Diffuses each vertex toward the mean of its 1-ring neighbors. With ``implicit_time_integration``
    the diffusion is solved implicitly (backward Euler, unconditionally stable); otherwise an
    explicit forward-Euler step is used (stable for ``lamb <= 1``). Mirrors
    [`trimesh.smoothing.filter_laplacian`][] but returns a new vertex array instead of mutating
    the mesh in place.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    lamb
        Diffusion speed. ``0`` leaves the mesh unchanged; larger values smooth faster.
    iterations
        Number of smoothing passes.
    implicit_time_integration
        If ``True`` solve ``((1 + lamb) I - lamb L) V' = V`` each pass via conjugate gradient.
        If ``False`` apply the explicit step ``V' = V + lamb (L V - V)``.
    volume_constraint
        If ``True`` rescale the mesh after each pass to preserve its initial volume, counteracting
        Laplacian shrinkage.
    laplacian_operator
        Optional precomputed row-stochastic operator (see
        [`laplacian`][triwarp.laplacian.laplacian]). Autogenerated (uniform weights) when ``None``.

    Returns
    -------
    wp.array[wp.vec3]
        Smoothed ``(n_vertices,)`` vertex positions on ``vertices.device``.

    Raises
    ------
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`filter_taubin`][triwarp.smoothing.filter_taubin]
    [`filter_humphrey`][triwarp.smoothing.filter_humphrey]
    [`filter_implicit_fairing`][triwarp.smoothing.filter_implicit_fairing]
    [`trimesh.smoothing.filter_laplacian`][]
    """
    require_same_device(vertices=vertices, faces=faces)
    # This six-line prologue is shared by five of the filters (here, ``filter_humphrey``,
    # ``filter_taubin``, ``filter_neighborhood_average`` and ``filter_mut_dif_laplacian``), but a
    # shared helper is not worth it: the two halves have no common consumer.
    # ``filter_implicit_fairing`` needs the guard and ``_as_vec3d`` but builds no operator,
    # ``filter_neighborhood_average`` needs the operator built ``symmetric=True``, and a helper
    # returning both would have to hand back an optional operator that three of the six callers
    # immediately unwrap.
    device = vertices.device
    n = int(vertices.shape[0])
    if n == 0 or iterations == 0:
        return wp.clone(vertices)

    operator = _resolved_operator(vertices, faces, laplacian_operator)
    positions = _as_vec3d(vertices)
    if volume_constraint:
        vol_ini = tw.measures.volume(positions, faces)
        # The *initial* centre of mass, computed once here rather than inside the loop -- see
        # _apply_volume_constraint's Notes for why the rescale has to stay anchored to this one
        # fixed point rather than the mesh's current (already-moved) centre of mass.
        _, center_f32, _ = tw.measures.moments(vertices, faces)
        center_ini = wp.vec3d(float(center_f32[0]), float(center_f32[1]), float(center_f32[2]))
    else:
        vol_ini = 0.0
        center_ini = wp.vec3d(0.0, 0.0, 0.0)

    if implicit_time_integration:
        system = _build_implicit_system(operator, lamb, n, device)
        components = _component_columns(n, device)
        solutions = _component_columns(n, device)
        # The operator is fixed for the whole flow -- only the right-hand side moves -- so the
        # batched solver state is built once outside the loop and re-reads both buffers on every
        # call, which is the shape its own docstring asks for.
        solver = twl.spd_column_solver(
            system, components, solutions, tol=twl.CG_TOLERANCE, maxiter=10 * n
        )
        component_rows = [components[column] for column in range(3)]
        for _ in range(iterations):
            wp.map(kernel_smoothing.extract_components, positions, out=component_rows)
            # Seed each column with its own right-hand side, which here *is* the current position
            # component. Only a volume-constrained pass makes the two differ, since the rescale
            # moves the positions after the previous solve wrote them.
            wp.copy(solutions, components)
            solver()
            wp.map(
                kernel_smoothing.combine_components,
                solutions[0],
                solutions[1],
                solutions[2],
                out=positions,
            )
            if volume_constraint:
                _apply_volume_constraint(positions, faces, vol_ini, center_ini)
    else:
        nxt = wp.empty(n, dtype=wp.vec3d, device=device)
        coeff = wp.float64(lamb)
        for _ in range(iterations):
            _diffuse_pass(operator, positions, coeff, nxt)
            positions, nxt = nxt, positions
            if volume_constraint:
                _apply_volume_constraint(positions, faces, vol_ini, center_ini)

    return _as_vec3(positions)

inflate 

inflate(
    vertices: array[vec3],
    faces: array[int32],
    pressure: float,
    *,
    iterations: int = 3,
    pre_smooth: bool = True,
    gradual: bool = True,
    lamb: float = 0.5,
) -> wp.array[wp.vec3]

Inflate a surface under uniform pressure: normal displacement alternating with smoothing.

The "balloon" flow. Each pass moves every vertex along its own vertex_normals at weighting="area" and then relaxes with one Laplacian pass, which is what keeps the triangles from shearing as the surface grows. Uses of it: puffing a thin shell out to a printable thickness, opening a collapsed scan, and supplying a starting surface a fitting loop can shrink back onto data.

pressure is an absolute distance per pass, in the mesh's own units, so scale it off something intrinsic -- a fraction of the mean edge length is the usual choice. With gradual pass k of n uses pressure * (k + 1) / n rather than the full amount, which is gentler on a mesh with fine triangles: the early passes let the relaxation redistribute before the later ones push hard.

Not the implicit formulation

The reference inflations displace and then solve an implicit Laplacian system, where this runs one explicit pass -- so the two produce different surfaces from the same pressure, and the tests compare the properties an inflation must have (volume grows, the displacement is normal-aligned) rather than positions. Pass the result through filter_implicit_fairing if you want that formulation; it costs a conjugate-gradient solve per pass, which is why it is not the default here.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
pressure float

Distance to move along the normal per pass. Negative deflates.

 required
iterations int

Number of displace-and-relax passes.

 3
pre_smooth bool

Run one relaxation pass before the first displacement. Worth it on a noisy input, where the normals are what the noise corrupts and displacing along them amplifies it.

 True
gradual bool

Ramp the pressure linearly across the passes instead of applying it in full each time.

 True
lamb float

Relaxation strength of each Laplacian pass, as filter_laplacian takes it.

 0.5

Returns:

Type Description
array[vec3]

Inflated positions on vertices.device. The connectivity is untouched, so faces remains valid.

Raises:

Type Description
ValueError

If iterations is negative.
RuntimeError

If vertices and faces are not all on one device.
See Also

offset_mesh The exact outward offset, through a signed distance field -- it changes the connectivity and cannot self-intersect, where this keeps the mesh and can. filter_laplacian The relaxation half of each pass. thicken_mesh Turns a surface into a solid shell, which is the other way to give it thickness.

Source code in triwarp/smoothing.py 
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
def inflate(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    pressure: float,
    *,
    iterations: int = 3,
    pre_smooth: bool = True,
    gradual: bool = True,
    lamb: float = 0.5,
) -> wp.array[wp.vec3]:
    """
    Inflate a surface under uniform pressure: normal displacement alternating with smoothing.

    The "balloon" flow. Each pass moves every vertex along its own
    [`vertex_normals`][triwarp.vertices.vertex_normals] at ``weighting="area"`` and then relaxes
    with one Laplacian pass, which is what keeps the triangles from shearing as the surface grows.
    Uses of it: puffing a thin shell out to a printable thickness, opening a collapsed scan, and
    supplying a starting surface a fitting loop can shrink back onto data.

    ``pressure`` is an absolute distance per pass, in the mesh's own units, so scale it off
    something intrinsic -- a fraction of the mean edge length is the usual choice. With ``gradual``
    pass ``k`` of ``n`` uses ``pressure * (k + 1) / n`` rather than the full amount, which is
    gentler on a mesh with fine triangles: the early passes let the relaxation redistribute
    before the later ones push hard.

    !!! note "Not the implicit formulation"
        The reference inflations displace and then solve an *implicit* Laplacian system, where
        this runs one explicit pass -- so the two produce different surfaces from the same
        ``pressure``, and the tests compare the properties an inflation must have (volume grows,
        the displacement is normal-aligned) rather than positions. Pass the result through
        [`filter_implicit_fairing`][triwarp.smoothing.filter_implicit_fairing] if you want that
        formulation; it costs a conjugate-gradient solve per pass, which is why it is not the
        default here.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    pressure
        Distance to move along the normal per pass. Negative deflates.
    iterations
        Number of displace-and-relax passes.
    pre_smooth
        Run one relaxation pass *before* the first displacement. Worth it on a noisy input, where
        the normals are what the noise corrupts and displacing along them amplifies it.
    gradual
        Ramp the pressure linearly across the passes instead of applying it in full each time.
    lamb
        Relaxation strength of each Laplacian pass, as
        [`filter_laplacian`][triwarp.smoothing.filter_laplacian] takes it.

    Returns
    -------
    wp.array[wp.vec3]
        Inflated positions on ``vertices.device``. The connectivity is untouched, so ``faces``
        remains valid.

    Raises
    ------
    ValueError
        If ``iterations`` is negative.
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`offset_mesh`][triwarp.levelset.offset_mesh]
        The *exact* outward offset, through a signed distance field -- it changes the connectivity
        and cannot self-intersect, where this keeps the mesh and can.
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian]
        The relaxation half of each pass.
    [`thicken_mesh`][triwarp.levelset.thicken_mesh]
        Turns a surface into a solid shell, which is the other way to give it thickness.
    """
    require_same_device(vertices=vertices, faces=faces)
    if iterations < 0:
        raise ValueError(f"iterations must be non-negative, got {iterations}")
    device = faces.device
    n_vertices = int(vertices.shape[0])
    if n_vertices == 0 or int(faces.shape[0]) == 0 or iterations == 0:
        return wp.clone(vertices)

    positions = wp.clone(vertices)
    # The relaxation operator is hoisted for the same reason the step kernel below is, and it is the
    # larger of the two: ``filter_laplacian`` builds the uniform operator per call, and at
    # ``equal_weight=True`` that operator is the mesh's topology, which no pass changes.
    operator = laplacian.laplacian(vertices, faces)
    if pre_smooth:
        positions = filter_laplacian(
            positions, faces, lamb, iterations=1, laplacian_operator=operator
        )
    # Hoisted once: the displacement is the same map every pass, so a per-iteration wrapper loop
    # should not re-derive it.
    step_kernel = wp.map(
        kernel_smoothing.step_along_normal,
        positions,
        positions,
        wp.float32(0.0),
        out=positions,
        return_kernel=True,
    )
    for step in range(iterations):
        amount = pressure * (step + 1) / iterations if gradual else pressure
        normals = tw.vertices.vertex_normals(positions, faces)
        displaced = wp.empty(n_vertices, dtype=wp.vec3, device=device)
        wp.launch(
            step_kernel,
            dim=n_vertices,
            inputs=[positions, normals, wp.float32(amount)],
            outputs=[displaced],
            device=device,
        )
        positions = filter_laplacian(
            displaced, faces, lamb, iterations=1, laplacian_operator=operator
        )
    return positions

filter_humphrey 

filter_humphrey(
    vertices: array[vec3],
    faces: array[int32],
    alpha: float = 0.1,
    beta: float = 0.5,
    iterations: int = 10,
    laplacian_operator: BsrMatrix[float32] | None = None,
) -> wp.array[wp.vec3]

Laplacian smoothing with Humphrey (HC) filtering (Vollmer et al.).

Each pass performs a Laplacian smoothing step and then pushes the vertices partway back toward their original and pre-smoothing positions, strongly reducing the shrinkage of plain Laplacian smoothing. Mirrors trimesh.smoothing.filter_humphrey.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
alpha float

Shrinkage control in [0, 1]: 0 ignores the original positions, 1 disables smoothing.

 0.1
beta float

Correction aggressiveness in [0, 1]: 0 disables smoothing, 1 is most aggressive.

 0.5
iterations int

Number of smoothing passes.

 10
laplacian_operator BsrMatrix[float32] | None

Optional precomputed row-stochastic operator (see laplacian). Autogenerated (uniform weights) when None.

 None

Returns:

Type Description
array[vec3]

Smoothed (n_vertices,) vertex positions on vertices.device.

Raises:

Type Description
RuntimeError

If vertices and faces are not all on one device.
See Also

filter_laplacian trimesh.smoothing.filter_humphrey

Source code in triwarp/smoothing.py 
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
def filter_humphrey(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    alpha: float = 0.1,
    beta: float = 0.5,
    iterations: int = 10,
    laplacian_operator: wps.BsrMatrix[wp.float32] | None = None,
) -> wp.array[wp.vec3]:
    """
    Laplacian smoothing with Humphrey (HC) filtering (Vollmer et al.).

    Each pass performs a Laplacian smoothing step and then pushes the vertices partway back
    toward their original and pre-smoothing positions, strongly reducing the shrinkage of plain
    Laplacian smoothing. Mirrors [`trimesh.smoothing.filter_humphrey`][].

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    alpha
        Shrinkage control in ``[0, 1]``: ``0`` ignores the original positions, ``1`` disables
        smoothing.
    beta
        Correction aggressiveness in ``[0, 1]``: ``0`` disables smoothing, ``1`` is most aggressive.
    iterations
        Number of smoothing passes.
    laplacian_operator
        Optional precomputed row-stochastic operator (see
        [`laplacian`][triwarp.laplacian.laplacian]). Autogenerated (uniform weights) when ``None``.

    Returns
    -------
    wp.array[wp.vec3]
        Smoothed ``(n_vertices,)`` vertex positions on ``vertices.device``.

    Raises
    ------
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian]
    [`trimesh.smoothing.filter_humphrey`][]
    """
    require_same_device(vertices=vertices, faces=faces)
    device = vertices.device
    n = int(vertices.shape[0])
    if n == 0 or iterations == 0:
        return wp.clone(vertices)

    operator = _resolved_operator(vertices, faces, laplacian_operator)
    positions = _as_vec3d(vertices)
    original = wp.clone(positions)

    lv = wp.empty(n, dtype=wp.vec3d, device=device)
    b = wp.empty(n, dtype=wp.vec3d, device=device)
    nxt = wp.empty(n, dtype=wp.vec3d, device=device)
    alpha64 = wp.float64(alpha)
    beta64 = wp.float64(beta)
    rows = (operator.offsets, operator.columns, operator.values)
    for _ in range(iterations):
        # Two passes rather than four: each applies the operator in the thread that consumes its
        # row. ``lv`` still crosses between them -- the update needs ``L.v`` *and* ``L.b`` -- but
        # ``L.b`` never leaves a register. ``positions`` doubles as the previous-iterate ``q`` (it
        # is only read this pass).
        wp.launch(
            kernel_smoothing.humphrey_residual_pass,
            dim=n,
            inputs=[*rows, positions, original, alpha64],
            outputs=[lv, b],
            device=device,
        )
        wp.launch(
            kernel_smoothing.humphrey_update_pass,
            dim=n,
            inputs=[*rows, lv, b, beta64],
            outputs=[nxt],
            device=device,
        )
        positions, nxt = nxt, positions

    return _as_vec3(positions)

filter_spikes 

filter_spikes(
    vertices: array[vec3],
    faces: array[int32],
    min_angle_sum: float,
    *,
    max_iter: int = 10,
    return_count: bool = False,
) -> wp.array[wp.vec3] | tuple[wp.array[wp.vec3], int]

Pull needle-like vertices back onto the surface, leaving every other vertex untouched.

A spike is a vertex whose incident corner angles sum to less than min_angle_sum. On a flat surface that sum is 2 * pi; the sharper the cone, the smaller it gets, and a vertex on a thin needle -- the classic scan and reconstruction artifact -- has a sum near zero. So the threshold is read in radians against a full turn, and a value like 0.5 * pi selects genuinely degenerate cones while leaving ordinary sharp features alone.

Each pass replaces the flagged vertices by the average of their closed 1-rings and nothing else moves: this is not a smoothing filter restricted to a region, it is a repair that touches only the vertices that fail the test. Flattening a spike changes its neighbours' angle sums, so the pass repeats -- a needle whose base is itself spiky needs more than one -- and stops as soon as a pass finds none.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
min_angle_sum float

Angle-sum threshold in radians. A vertex is a spike when its incident angles sum below it. Compare against 2 * pi, the flat value.

 required
max_iter int

Cap on the number of passes.

 10
return_count bool

If True, also return flattened.

 False

Returns:

Name Type Description
vertices array[vec3]

Repaired positions on vertices.device. Connectivity is untouched, so faces stays valid.
flattened (int, optional)

Present when return_count=True. How many vertex moves were made, summed over the passes -- so a vertex fixed twice counts twice, which makes it a report rather than a count of anything in the output. Zero means nothing was flagged and the buffer is the input's. The pass loop stops itself as soon as a pass finds no spike, so the default return is a bare position buffer that drops into a filter chain like every other function here.

Raises:

Type Description
ValueError

If max_iter is negative.
RuntimeError

If vertices and faces are not all on one device.
See Also

filter_neighborhood_average The unrestricted version of the move this makes, over every vertex. vertex_defects 2 * pi minus the same angle sum, which is what the threshold is read against. validation.face_defective_mask Flags the faces a spike produces, where this flags the vertex itself.

Source code in triwarp/smoothing.py 
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
def filter_spikes(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    min_angle_sum: float,
    *,
    max_iter: int = 10,
    return_count: bool = False,
) -> wp.array[wp.vec3] | tuple[wp.array[wp.vec3], int]:
    """
    Pull needle-like vertices back onto the surface, leaving every other vertex untouched.

    A *spike* is a vertex whose incident corner angles sum to less than ``min_angle_sum``. On a
    flat surface that sum is ``2 * pi``; the sharper the cone, the smaller it gets, and a vertex on
    a thin needle -- the classic scan and reconstruction artifact -- has a sum near zero. So the
    threshold is read in radians against a full turn, and a value like ``0.5 * pi`` selects
    genuinely degenerate cones while leaving ordinary sharp features alone.

    Each pass replaces the flagged vertices by the average of their closed 1-rings and **nothing
    else moves**: this is not a smoothing filter restricted to a region, it is a repair that touches
    only the vertices that fail the test. Flattening a spike changes its neighbours' angle sums, so
    the pass repeats -- a needle whose base is itself spiky needs more than one -- and stops as soon
    as a pass finds none.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    min_angle_sum
        Angle-sum threshold in radians. A vertex is a spike when its incident angles sum below it.
        Compare against ``2 * pi``, the flat value.
    max_iter
        Cap on the number of passes.
    return_count
        If ``True``, also return ``flattened``.

    Returns
    -------
    vertices : wp.array[wp.vec3]
        Repaired positions on ``vertices.device``. Connectivity is untouched, so ``faces`` stays
        valid.
    flattened : int, optional
        Present when ``return_count=True``. How many vertex moves were made, summed over the passes
        -- so a vertex fixed twice counts twice, which makes it a report rather than a count of
        anything in the output. Zero means nothing was flagged and the buffer is the input's. The
        pass loop stops itself as soon as a pass finds no spike, so the default return is a bare
        position buffer that drops into a filter chain like every other function here.

    Raises
    ------
    ValueError
        If ``max_iter`` is negative.
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`filter_neighborhood_average`][triwarp.smoothing.filter_neighborhood_average]
        The unrestricted version of the move this makes, over every vertex.
    [`vertex_defects`][triwarp.vertices.vertex_defects]
        ``2 * pi`` minus the same angle sum, which is what the threshold is read against.
    [`validation.face_defective_mask`][triwarp.validation.face_defective_mask]
        Flags the *faces* a spike produces, where this flags the vertex itself.
    """
    require_same_device(vertices=vertices, faces=faces)
    if max_iter < 0:
        raise ValueError(f"max_iter must be non-negative, got {max_iter}")
    device = faces.device
    n_vertices = int(vertices.shape[0])
    if n_vertices == 0 or int(faces.shape[0]) == 0 or max_iter == 0:
        cloned = wp.clone(vertices)
        return (cloned, 0) if return_count else cloned

    positions = wp.clone(vertices)
    flattened = 0
    spikes = wp.empty(n_vertices, dtype=wp.bool, device=device)
    # Built once and handed to every pass. ``filter_neighborhood_average`` would build it per call,
    # and at ``equal_weight=True`` -- its default, and what ``_resolved_operator`` asks for -- the
    # operator reads no positions at all: every off-diagonal weight is ``1`` before the row
    # normalization, so it is the mesh's *topology*, which no pass changes. Hoisting it is therefore
    # exactly equivalent, and the cost of a build is flat in the mesh size.
    operator = laplacian.laplacian(vertices, faces, symmetric=True)
    for _ in range(max_iter):
        defects = tw.vertices.vertex_defects(
            n_vertices, faces, tw.triangles.face_angles(positions, faces)
        )
        # The spike test is ``angle_sum < min_angle_sum``, applied to the angle *defect*
        # ``2 * pi - angle_sum`` instead: the condition becomes ``defect > 2 * pi - min_angle_sum``,
        # and the defect is the quantity ``vertices.vertex_defects`` already returns -- re-deriving
        # the sum would mean scattering the same corner angles a second time.
        wp.map(kernel_array.greater, defects, wp.float32(2.0 * math.pi - min_angle_sum), out=spikes)
        # One readback per pass, and it is the stopping test: whether any vertex is still a spike is
        # a device-side fact that a Python loop cannot branch on otherwise.
        n_spikes = int(tw.reduce.sum(tw.array.astype(spikes, wp.int32)))
        if n_spikes == 0:
            break
        smoothed = filter_neighborhood_average(
            positions, faces, iterations=1, laplacian_operator=operator
        )
        wp.map(kernel_smoothing.select_position, smoothed, positions, spikes, out=positions)
        flattened += n_spikes
    return (positions, flattened) if return_count else positions

equalize_triangle_areas 

equalize_triangle_areas(
    vertices: array[vec3],
    faces: array[int32],
    iterations: int = 1,
    force: float = 0.5,
    no_shrinkage: bool = False,
    region: array[bool] | None = None,
    max_displacement: float | None = None,
    *,
    vertex_faces: tuple[array[int32], array[int32]]
    | None = None,
) -> wp.array[wp.vec3]

Even out triangle areas by moving vertices, without touching the connectivity.

Every other filter in this module drives the surface toward a Laplacian residual of zero, which equalizes edge directions and says nothing about area: a patch of long thin triangles beside a patch of fat ones is already Laplacian-smooth. This one minimizes the summed squared areas of each vertex's incident triangles instead, which is the objective that actually redistributes them -- a vertex sitting close to one of its opposite edges is pushed away from it.

Twice a triangle's area is |(x - p) x (q - p)| for the free vertex x over its opposite edge (p, q), so the per-vertex objective is a sum of quadratic forms and its minimum is one 3x3 solve, taken in float64. Each pass steps a fraction force of the way there, from the previous pass's positions, so the whole mesh moves at once and the result does not depend on a vertex ordering.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
iterations int

Number of passes. 0 returns a copy.

 1
force float

Fraction of the way to the per-vertex minimum each pass moves. 1.0 jumps straight there and oscillates on an irregular mesh; the default is the usual under-relaxation.

 0.5
no_shrinkage bool

Constrain each vertex to the tangent plane through its current position, so it slides across the surface instead of sinking into it. The unconstrained minimum of squared area pulls the whole 1-ring inward, so leave this on when the shape matters and off when only the triangle quality does.

 False
region array[bool] | None

(n_vertices,) boolean mask of the vertices allowed to move. None moves every one. Vertices outside it stay exactly where they are, and are still read as neighbours.

 None
max_displacement float | None

Clamp every vertex into a ball of this radius around its input position, applied after each pass. None leaves the relaxation unbounded. Use it when the surface has to stay within a tolerance of the scan it came from.

 None
vertex_faces tuple[array[int32], array[int32]] | None

Optional precomputed vertex_face_adjacency as (vertex_faces, offsets). Depends on the connectivity alone -- which this filter never changes -- so one CSR serves every call over the same mesh, and Trimesh.vertex_face_adjacency has it cached.

 None

Returns:

Type Description
array[vec3]

Relaxed (n_vertices,) positions on vertices.device. Connectivity is untouched, so faces stays valid.

Raises:

Type Description
ValueError

If iterations is negative, if max_displacement is negative, or if region is not a length-n_vertices wp.bool array.
RuntimeError

If vertices, faces, region and vertex_faces are not all on one device.
See Also

relax_keep_volume The other member of this group: a uniform relax that does not shrink the shape. filter_laplacian The plain diffusion this is the area-driven alternative to. isotropic_remesh Equalizes edge lengths, by changing the connectivity as well.

Notes

It equalizes areas, not shapes

Nothing here bounds a triangle's aspect ratio, and a long thin triangle can have exactly the right area. For shape, flip the triangulation (flip_by_objective) or remesh (isotropic_remesh); this moves vertices only.

Source code in triwarp/smoothing.py 
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
def equalize_triangle_areas(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    iterations: int = 1,
    force: float = 0.5,
    no_shrinkage: bool = False,
    region: wp.array[wp.bool] | None = None,
    max_displacement: float | None = None,
    *,
    vertex_faces: tuple[wp.array[wp.int32], wp.array[wp.int32]] | None = None,
) -> wp.array[wp.vec3]:
    """
    Even out triangle *areas* by moving vertices, without touching the connectivity.

    Every other filter in this module drives the surface toward a Laplacian residual of zero, which
    equalizes *edge* directions and says nothing about area: a patch of long thin triangles beside a
    patch of fat ones is already Laplacian-smooth. This one minimizes the summed squared areas of
    each vertex's incident triangles instead, which is the objective that actually redistributes
    them -- a vertex sitting close to one of its opposite edges is pushed away from it.

    Twice a triangle's area is ``|(x - p) x (q - p)|`` for the free vertex ``x`` over its opposite
    edge ``(p, q)``, so the per-vertex objective is a sum of quadratic forms and its minimum is one
    3x3 solve, taken in ``float64``. Each pass steps a fraction ``force`` of the way there, from the
    *previous* pass's positions, so the whole mesh moves at once and the result does not depend on
    a vertex ordering.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    iterations
        Number of passes. ``0`` returns a copy.
    force
        Fraction of the way to the per-vertex minimum each pass moves. ``1.0`` jumps straight there
        and oscillates on an irregular mesh; the default is the usual under-relaxation.
    no_shrinkage
        Constrain each vertex to the tangent plane through its current position, so it slides
        across the surface instead of sinking into it. The unconstrained minimum of *squared* area
        pulls the whole 1-ring inward, so leave this on when the shape matters and off when only
        the triangle quality does.
    region
        ``(n_vertices,)`` boolean mask of the vertices allowed to move. ``None`` moves every one.
        Vertices outside it stay exactly where they are, and are still read as neighbours.
    max_displacement
        Clamp every vertex into a ball of this radius around its **input** position, applied after
        each pass. ``None`` leaves the relaxation unbounded. Use it when the surface has to stay
        within a tolerance of the scan it came from.
    vertex_faces
        Optional precomputed [`vertex_face_adjacency`][triwarp.adjacency.vertex_face_adjacency] as
        ``(vertex_faces, offsets)``. Depends on the connectivity alone -- which this filter never
        changes -- so one CSR serves every call over the same mesh, and
        [`Trimesh.vertex_face_adjacency`][triwarp.mesh.Trimesh.vertex_face_adjacency] has it cached.

    Returns
    -------
    wp.array[wp.vec3]
        Relaxed ``(n_vertices,)`` positions on ``vertices.device``. Connectivity is untouched, so
        ``faces`` stays valid.

    Raises
    ------
    ValueError
        If ``iterations`` is negative, if ``max_displacement`` is negative, or if ``region`` is not
        a length-``n_vertices`` ``wp.bool`` array.
    RuntimeError
        If ``vertices``, ``faces``, ``region`` and ``vertex_faces`` are not all on one device.

    See Also
    --------
    [`relax_keep_volume`][triwarp.smoothing.relax_keep_volume]
        The other member of this group: a uniform relax that does not shrink the shape.
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian]
        The plain diffusion this is the area-driven alternative to.
    [`isotropic_remesh`][triwarp.remesh.isotropic_remesh]
        Equalizes edge *lengths*, by changing the connectivity as well.

    Notes
    -----
    !!! note "It equalizes areas, not shapes"
        Nothing here bounds a triangle's aspect ratio, and a long thin triangle can have exactly the
        right area. For shape, flip the triangulation
        ([`flip_by_objective`][triwarp.remesh.flip_by_objective]) or remesh
        ([`isotropic_remesh`][triwarp.remesh.isotropic_remesh]); this moves vertices only.
    """
    require_same_device(vertices=vertices, faces=faces, region=region, vertex_faces=vertex_faces)
    device = vertices.device
    n_vertices = int(vertices.shape[0])
    flags, limit = _relaxation_state(vertices, iterations, region, max_displacement)
    if flags is None:
        return wp.clone(vertices)

    vf_indices, offsets = (
        vertex_faces
        if vertex_faces is not None
        else tw.adjacency.vertex_face_adjacency(faces, n_vertices=n_vertices)
    )
    positions = wp.clone(vertices)
    normals = wp.zeros(n_vertices, dtype=wp.vec3, device=device)
    nxt = wp.empty(n_vertices, dtype=wp.vec3, device=device)
    for _ in range(iterations):
        if no_shrinkage:
            # Recomputed per pass rather than hoisted: the tangent plane the solve is confined to is
            # the one through the vertex *now*, and after a pass that is a different plane.
            normals = tw.vertices.vertex_normals(positions, faces)
        wp.launch(
            kernel_smoothing.equalize_area_step,
            dim=n_vertices,
            inputs=[
                positions,
                faces,
                offsets,
                vf_indices,
                normals,
                flags,
                vertices,
                wp.float32(force),
                no_shrinkage,
                limit,
                nxt,
            ],
            device=device,
        )
        positions, nxt = nxt, positions
    return positions

relax_keep_volume 

relax_keep_volume(
    vertices: array[vec3],
    faces: array[int32],
    iterations: int = 1,
    force: float = 0.5,
    region: array[bool] | None = None,
    max_displacement: float | None = None,
) -> wp.array[wp.vec3]

Uniform relaxation with the shrinkage taken out locally, rather than rescaled away.

A plain 1-ring average moves every vertex toward the inside of the surface it curves around, so a closed shape loses volume with every pass. filter_laplacian's volume_constraint answers that by rescaling the whole mesh at the end, which restores the total and is wrong everywhere the shrinkage was not uniform. This answers it per vertex instead: the relax displacement is computed as a field, and each vertex then has its own neighbourhood's average displacement subtracted from its own.

That difference is the whole method. A translation shared by a neighbourhood -- which is what shrinkage locally is -- cancels exactly, while the high-frequency part of the displacement, which is the noise the relax was for, survives untouched. No global rescale, and a mesh with a flat region and a curved one is treated correctly in both.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
iterations int

Number of passes. 0 returns a copy.

 1
force float

Fraction of the way to the 1-ring average each pass moves before the correction.

 0.5
region array[bool] | None

(n_vertices,) boolean mask of the vertices allowed to move. None moves every one. A neighbour outside the region contributes its position to the average but no displacement to the correction, so a vertex on the region's edge is corrected by less than an interior one -- which is what stops the region from tearing away from the rest.

 None
max_displacement float | None

Clamp every vertex into a ball of this radius around its input position, applied after each pass. None leaves the relaxation unbounded.

 None

Returns:

Type Description
array[vec3]

Relaxed (n_vertices,) positions on vertices.device.

Raises:

Type Description
ValueError

If iterations is negative, if max_displacement is negative, or if region is not a length-n_vertices wp.bool array.
RuntimeError

If vertices, faces and region are not all on one device.
See Also

filter_laplacian The plain diffusion, whose volume_constraint is the global rescale this replaces. filter_taubin The other anti-shrinkage answer: alternate a shrinking pass with an inflating one. equalize_triangle_areas

Source code in triwarp/smoothing.py 
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
def relax_keep_volume(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    iterations: int = 1,
    force: float = 0.5,
    region: wp.array[wp.bool] | None = None,
    max_displacement: float | None = None,
) -> wp.array[wp.vec3]:
    """
    Uniform relaxation with the shrinkage taken out **locally**, rather than rescaled away.

    A plain 1-ring average moves every vertex toward the inside of the surface it curves around,
    so a closed shape loses volume with every pass.
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian]'s ``volume_constraint`` answers that
    by rescaling the whole mesh at the end, which restores the total and is wrong everywhere the
    shrinkage was not uniform. This answers it per vertex instead: the relax displacement is
    computed as a *field*, and each vertex then has its own neighbourhood's **average
    displacement** subtracted from its own.

    That difference is the whole method. A translation shared by a neighbourhood -- which is what
    shrinkage locally is -- cancels exactly, while the high-frequency part of the displacement,
    which is the noise the relax was for, survives untouched. No global rescale, and a mesh with a
    flat region and a curved one is treated correctly in both.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    iterations
        Number of passes. ``0`` returns a copy.
    force
        Fraction of the way to the 1-ring average each pass moves before the correction.
    region
        ``(n_vertices,)`` boolean mask of the vertices allowed to move. ``None`` moves every one.
        A neighbour outside the region contributes its position to the average but no displacement
        to the correction, so a vertex on the region's edge is corrected by less than an interior
        one -- which is what stops the region from tearing away from the rest.
    max_displacement
        Clamp every vertex into a ball of this radius around its **input** position, applied after
        each pass. ``None`` leaves the relaxation unbounded.

    Returns
    -------
    wp.array[wp.vec3]
        Relaxed ``(n_vertices,)`` positions on ``vertices.device``.

    Raises
    ------
    ValueError
        If ``iterations`` is negative, if ``max_displacement`` is negative, or if ``region`` is not
        a length-``n_vertices`` ``wp.bool`` array.
    RuntimeError
        If ``vertices``, ``faces`` and ``region`` are not all on one device.

    See Also
    --------
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian]
        The plain diffusion, whose ``volume_constraint`` is the *global* rescale this replaces.
    [`filter_taubin`][triwarp.smoothing.filter_taubin]
        The other anti-shrinkage answer: alternate a shrinking pass with an inflating one.
    [`equalize_triangle_areas`][triwarp.smoothing.equalize_triangle_areas]
    """
    require_same_device(vertices=vertices, faces=faces, region=region)
    device = vertices.device
    n_vertices = int(vertices.shape[0])
    flags, limit = _relaxation_state(vertices, iterations, region, max_displacement)
    if flags is None:
        return wp.clone(vertices)

    unique_edges, _ = tw.edges.edges_unique(faces, n_vertices=n_vertices)
    adjacency = tw.graph.edges_to_csr(n_vertices, unique_edges)
    offsets, columns = adjacency.offsets, adjacency.columns
    positions = wp.clone(vertices)
    push = wp.empty(n_vertices, dtype=wp.vec3, device=device)
    nxt = wp.empty(n_vertices, dtype=wp.vec3, device=device)
    for _ in range(iterations):
        wp.launch(
            kernel_smoothing.ring_push_forces,
            dim=n_vertices,
            inputs=[positions, offsets, columns, flags, wp.float32(force), push],
            device=device,
        )
        wp.launch(
            kernel_smoothing.apply_push_keeping_volume,
            dim=n_vertices,
            inputs=[positions, offsets, columns, flags, push, vertices, limit, nxt],
            device=device,
        )
        positions, nxt = nxt, positions
    return positions

relax_approx 

relax_approx(
    vertices: array[vec3],
    faces: array[int32],
    dilate_radius: float,
    iterations: int = 1,
    force: float = 0.5,
    fit: Literal["planar", "quadric"] = "planar",
    region: array[bool] | None = None,
    max_displacement: float | None = None,
) -> wp.array[wp.vec3]

Move each vertex onto a surface fitted to its neighbourhood, rather than toward its 1-ring.

Every diffusion filter here is a 1-ring operator, so the largest feature it can see is one edge across and its idea of "where the surface is" is the average of six neighbours. This fits a local surface -- a least-squares plane, or a quadric graph over that plane -- to every vertex inside a geodesic ball of dilate_radius and moves the vertex onto it. The ball is what makes it different in kind: noise smaller than the radius is averaged out in one pass instead of being diffused away over many, and the fitted surface is not pulled off the shape by the vertex being fitted.

fit="planar" flattens: it is the strongest of the two and will take the curvature out of a genuinely curved surface if the radius is large. fit="quadric" keeps curvature, because a curved surface is in its model space, and is the one to use on anything but a plate.

The neighbourhood is a geodesic ball (geodesic_ball), not a Euclidean one, so the opposite wall of a thin tube is excluded even when it is close -- which is exactly the case that corrupts a fit.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
dilate_radius float

Radius of the neighbourhood ball, in world units. It has no scale-free default: below one edge length the ball is the vertex alone and the call is a no-op, and the useful range starts at a small multiple of the mean edge length (edges_unique_length).

 required
iterations int

Number of passes. 0 returns a copy. The neighbourhoods are built once, from the input positions, and reused -- the ball's membership is a topological choice and rebuilding it per pass would cost more than the fit.

 1
force float

Fraction of the way to the fitted surface each pass moves.

 0.5
fit Literal['planar', 'quadric']

"planar" for a least-squares plane, "quadric" for a 6-coefficient quadric graph over it.

 'planar'
region array[bool] | None

(n_vertices,) boolean mask of the vertices allowed to move. None moves every one. Vertices outside it still populate their neighbours' balls.

 None
max_displacement float | None

Clamp every vertex into a ball of this radius around its input position, applied after each pass. None leaves the relaxation unbounded.

 None

Returns:

Type Description
array[vec3]

Relaxed (n_vertices,) positions on vertices.device.

Raises:

Type Description
ValueError

If iterations is negative, if dilate_radius is not positive, if max_displacement is negative, if fit is not one of the two names, or if region is not a length-n_vertices wp.bool array.
RuntimeError

If vertices, faces and region are not all on one device.
See Also

filter_laplacian The 1-ring diffusion this replaces with a neighbourhood fit. principal_curvature Fits a quadric over the same geodesic balls, to measure rather than to move. geodesic_ball

Notes

A vertex with too small a ball does not move

A plane fit needs three independent points and the quadric six, so a vertex whose ball holds fewer than six vertices is left exactly where it is rather than fitted to whatever it has. On a mesh whose edges are longer than dilate_radius that is every vertex and the call returns the input -- check the result moved before concluding the parameters were right.

Source code in triwarp/smoothing.py 
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
def relax_approx(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    dilate_radius: float,
    iterations: int = 1,
    force: float = 0.5,
    fit: Literal["planar", "quadric"] = "planar",
    region: wp.array[wp.bool] | None = None,
    max_displacement: float | None = None,
) -> wp.array[wp.vec3]:
    """
    Move each vertex onto a surface fitted to its **neighbourhood**, rather than toward its 1-ring.

    Every diffusion filter here is a 1-ring operator, so the largest feature it can see is one edge
    across and its idea of "where the surface is" is the average of six neighbours. This fits a
    local surface -- a least-squares plane, or a quadric graph over that plane -- to every vertex
    inside a geodesic ball of ``dilate_radius`` and moves the vertex onto it. The ball is what makes
    it different in kind: noise smaller than the radius is averaged out in one pass instead of being
    diffused away over many, and the fitted surface is not pulled off the shape by the vertex being
    fitted.

    ``fit="planar"`` flattens: it is the strongest of the two and will take the curvature out of a
    genuinely curved surface if the radius is large. ``fit="quadric"`` keeps curvature, because a
    curved surface *is* in its model space, and is the one to use on anything but a plate.

    The neighbourhood is a **geodesic** ball ([`geodesic_ball`][triwarp.neighbors.geodesic_ball]),
    not a Euclidean one, so the opposite wall of a thin tube is excluded even when it is close --
    which is exactly the case that corrupts a fit.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    dilate_radius
        Radius of the neighbourhood ball, in world units. It has no scale-free default: below one
        edge length the ball is the vertex alone and the call is a no-op, and the useful range
        starts at a small multiple of the mean edge length
        ([`edges_unique_length`][triwarp.edges.edges_unique_length]).
    iterations
        Number of passes. ``0`` returns a copy. The neighbourhoods are built **once**, from the
        input positions, and reused -- the ball's membership is a topological choice and rebuilding
        it per pass would cost more than the fit.
    force
        Fraction of the way to the fitted surface each pass moves.
    fit
        ``"planar"`` for a least-squares plane, ``"quadric"`` for a 6-coefficient quadric graph over
        it.
    region
        ``(n_vertices,)`` boolean mask of the vertices allowed to move. ``None`` moves every one.
        Vertices outside it still populate their neighbours' balls.
    max_displacement
        Clamp every vertex into a ball of this radius around its **input** position, applied after
        each pass. ``None`` leaves the relaxation unbounded.

    Returns
    -------
    wp.array[wp.vec3]
        Relaxed ``(n_vertices,)`` positions on ``vertices.device``.

    Raises
    ------
    ValueError
        If ``iterations`` is negative, if ``dilate_radius`` is not positive, if
        ``max_displacement`` is negative, if ``fit`` is not one of the two names, or if ``region``
        is not a length-``n_vertices`` ``wp.bool`` array.
    RuntimeError
        If ``vertices``, ``faces`` and ``region`` are not all on one device.

    See Also
    --------
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian]
        The 1-ring diffusion this replaces with a neighbourhood fit.
    [`principal_curvature`][triwarp.curvature.principal_curvature]
        Fits a quadric over the same geodesic balls, to measure rather than to move.
    [`geodesic_ball`][triwarp.neighbors.geodesic_ball]

    Notes
    -----
    !!! note "A vertex with too small a ball does not move"
        A plane fit needs three independent points and the quadric six, so a vertex whose ball holds
        fewer than six vertices is left exactly where it is rather than fitted to whatever it has.
        On a mesh whose edges are longer than ``dilate_radius`` that is *every* vertex and the call
        returns the input -- check the result moved before concluding the parameters were right.
    """
    require_same_device(vertices=vertices, faces=faces, region=region)
    if fit not in ("planar", "quadric"):
        raise ValueError(f"fit must be 'planar' or 'quadric', got {fit!r}")
    if dilate_radius <= 0.0:
        raise ValueError(f"dilate_radius must be positive, got {dilate_radius}")
    device = vertices.device
    n_vertices = int(vertices.shape[0])
    flags, limit = _relaxation_state(vertices, iterations, region, max_displacement)
    if flags is None:
        return wp.clone(vertices)

    # Built once from the input positions: which vertices are in the ball is a decision about the
    # surface's connectivity, and the relaxation moves everything by less than the radius anyway.
    # ``min_count=1`` disables the nearest-neighbour backfill: a ball too small to fit is a
    # documented no-op here, and backfilling would quietly fit something else instead.
    neighbor_indices, neighbor_offsets, _ = tw.neighbors.geodesic_ball(
        vertices, faces, dilate_radius, min_count=1
    )
    positions = wp.clone(vertices)
    nxt = wp.empty(n_vertices, dtype=wp.vec3, device=device)
    for _ in range(iterations):
        wp.launch(
            kernel_smoothing.relax_approx_step,
            dim=n_vertices,
            inputs=[
                positions,
                neighbor_indices,
                neighbor_offsets,
                flags,
                vertices,
                wp.float32(force),
                fit == "quadric",
                limit,
                nxt,
            ],
            device=device,
        )
        positions, nxt = nxt, positions
    return positions

filter_taubin 

filter_taubin(
    vertices: array[vec3],
    faces: array[int32],
    lamb: float = 0.5,
    nu: float = 0.5,
    iterations: int = 10,
    laplacian_operator: BsrMatrix[float32] | None = None,
    *,
    recompute: bool = False,
) -> wp.array[wp.vec3]

Taubin lambda/nu mesh smoothing (Vollmer et al.).

Alternates a shrinking diffusion step (+lamb) on even passes with an inflating step (-nu) on odd passes, so the low-frequency shape is preserved while high-frequency noise is removed with little volume loss. Mirrors trimesh.smoothing.filter_taubin.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
lamb float

Shrinking step size in [0, 1].

 0.5
nu float

Inflating step size in [0, 1]; typically chosen so 0 < 1/lamb - 1/nu < 0.1.

 0.5
iterations int

Number of smoothing passes (shrink and inflate alternate by pass index).

 10
laplacian_operator BsrMatrix[float32] | None

Optional precomputed row-stochastic operator (see laplacian). Autogenerated (uniform weights) when None. Passing one and setting recompute is a contradiction and raises.

 None
recompute bool

Reassemble the inverse-distance operator from the current positions before every pass, instead of applying one operator throughout. Off by default; see the Notes — it is a different filter and markedly more expensive per pass.

 False

Returns:

Type Description
array[vec3]

Smoothed (n_vertices,) vertex positions on vertices.device.

Raises:

Type Description
ValueError

If recompute is set together with laplacian_operator, which would be ignored.
RuntimeError

If vertices and faces are not all on one device.
See Also

filter_laplacian trimesh.smoothing.filter_taubin laplacian The operator this applies, and the equal_weight switch recompute pins to False.

Notes

recompute is a different filter, not a tuning of this one: pytorch3d.ops.taubin_smoothing rebuilds its inverse-distance operator from the current geometry before each half-pass, and matching that convention is what makes recompute=True agree with it, where the fixed operator does not. The default is unchanged and stays the fixed operator, which is trimesh's filter and this group's oracle.

It is pinned to the inverse-distance weighting rather than taking a weighting of its own, and that is not a simplification: the uniform operator is 1 / degree off the diagonal, a function of the connectivity alone, which a position filter never changes -- so recomputing it would rebuild the identical matrix every pass and cost without doing anything. Only the geometry-dependent branch has anything to recompute.

recompute reassembles a sparse operator every pass instead of once, which is markedly more expensive; the connectivity never changes, though, so laplacian's edges keyword lets every pass share one edges_unique call rather than re-deriving it each time.

Source code in triwarp/smoothing.py 
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
def filter_taubin(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    lamb: float = 0.5,
    nu: float = 0.5,
    iterations: int = 10,
    laplacian_operator: wps.BsrMatrix[wp.float32] | None = None,
    *,
    recompute: bool = False,
) -> wp.array[wp.vec3]:
    """
    Taubin lambda/nu mesh smoothing (Vollmer et al.).

    Alternates a shrinking diffusion step (``+lamb``) on even passes with an inflating step
    (``-nu``) on odd passes, so the low-frequency shape is preserved while high-frequency noise is
    removed with little volume loss. Mirrors [`trimesh.smoothing.filter_taubin`][].

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    lamb
        Shrinking step size in ``[0, 1]``.
    nu
        Inflating step size in ``[0, 1]``; typically chosen so ``0 < 1/lamb - 1/nu < 0.1``.
    iterations
        Number of smoothing passes (shrink and inflate alternate by pass index).
    laplacian_operator
        Optional precomputed row-stochastic operator (see
        [`laplacian`][triwarp.laplacian.laplacian]). Autogenerated (uniform weights) when ``None``.
        Passing one *and* setting ``recompute`` is a contradiction and raises.
    recompute
        Reassemble the inverse-distance operator from the *current* positions before every pass,
        instead of applying one operator throughout. Off by default; see the ``Notes`` — it is a
        different filter and markedly more expensive per pass.

    Returns
    -------
    wp.array[wp.vec3]
        Smoothed ``(n_vertices,)`` vertex positions on ``vertices.device``.

    Raises
    ------
    ValueError
        If ``recompute`` is set together with ``laplacian_operator``, which would be ignored.
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian]
    [`trimesh.smoothing.filter_taubin`][]
    [`laplacian`][triwarp.laplacian.laplacian]
        The operator this applies, and the ``equal_weight`` switch ``recompute`` pins to ``False``.

    Notes
    -----
    ``recompute`` is a different filter, not a tuning of this one:
    ``pytorch3d.ops.taubin_smoothing`` rebuilds its inverse-distance operator from the current
    geometry before each half-pass, and matching that convention is what makes ``recompute=True``
    agree with it, where the fixed operator does not. The default is unchanged and stays the fixed
    operator, which is trimesh's filter and this group's oracle.

    It is pinned to the **inverse-distance** weighting rather than taking a weighting of its own,
    and that is not a simplification: the uniform operator is ``1 / degree`` off the diagonal, a
    function of the connectivity alone, which a position filter never changes -- so recomputing it
    would rebuild the identical matrix every pass and cost without doing anything. Only the
    geometry-dependent branch has anything to recompute.

    ``recompute`` reassembles a sparse operator every pass instead of once, which is markedly more
    expensive; the connectivity never changes, though, so ``laplacian``'s ``edges`` keyword lets
    every pass share one ``edges_unique`` call rather than re-deriving it each time.
    """
    require_same_device(vertices=vertices, faces=faces)
    if recompute and laplacian_operator is not None:
        raise ValueError("recompute rebuilds the operator each pass; do not also pass one")
    device = vertices.device
    n = int(vertices.shape[0])
    if n == 0 or iterations == 0:
        return wp.clone(vertices)

    operator = None if recompute else _resolved_operator(vertices, faces, laplacian_operator)
    # The recompute path reassembles the operator every pass from *moved* positions, but over
    # connectivity that never changes -- so the unique-edge set is derived once here rather than
    # inside the loop, which would otherwise repay the same derivation once per iteration for an
    # identical answer.
    recompute_edges = tw.edges.edges_unique(faces, n_vertices=n)[0] if recompute else None
    positions = _as_vec3d(vertices)
    nxt = wp.empty(n, dtype=wp.vec3d, device=device)
    for index in range(iterations):
        pass_operator = (
            laplacian.laplacian(
                _as_vec3(positions),
                faces,
                equal_weight=False,
                edges=recompute_edges,
                # ``recompute_edges`` was derived from this same ``faces``/``n`` a few lines up and
                # never changes across the loop, so re-checking it every pass would only re-pay the
                # per-iteration cost this hoist exists to remove.
                validate=False,
            )
            if operator is None
            else operator
        )
        coeff = lamb if index % 2 == 0 else -nu
        _diffuse_pass(pass_operator, positions, wp.float64(coeff), nxt)
        positions, nxt = nxt, positions

    return _as_vec3(positions)

filter_neighborhood_average 

filter_neighborhood_average(
    vertices: array[vec3],
    faces: array[int32],
    iterations: int = 10,
    laplacian_operator: BsrMatrix[float32] | None = None,
) -> wp.array[wp.vec3]

Uniform neighborhood-average smoothing over the closed 1-ring.

Each pass replaces every vertex by the equal-weight average of itself and its direct neighbors, v' = (v + sum(v_neighbors)) / (1 + degree). Unlike filter_laplacian (which averages the open 1-ring, excluding the vertex, and scales the result by lamb) this is the simplest possible smoother, with no diffusion coefficient. It is the CPU/GPU-portable analogue of Open3D's TriangleMesh.filter_smooth_simple.

Isolated vertices (degree 0) are left unchanged.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
iterations int

Number of smoothing passes.

 10
laplacian_operator BsrMatrix[float32] | None

Optional precomputed row-stochastic operator (see laplacian). Autogenerated (uniform weights) when None.

 None

Returns:

Type Description
array[vec3]

Smoothed (n_vertices,) vertex positions on vertices.device.

Raises:

Type Description
RuntimeError

If vertices and faces are not all on one device.
See Also

filter_laplacian filter_taubin

Source code in triwarp/smoothing.py 
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
def filter_neighborhood_average(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    iterations: int = 10,
    laplacian_operator: wps.BsrMatrix[wp.float32] | None = None,
) -> wp.array[wp.vec3]:
    """
    Uniform neighborhood-average smoothing over the closed 1-ring.

    Each pass replaces every vertex by the equal-weight average of itself and its direct
    neighbors, ``v' = (v + sum(v_neighbors)) / (1 + degree)``. Unlike
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian] (which averages the open 1-ring,
    excluding the vertex, and scales the result by ``lamb``) this is the simplest possible
    smoother, with no diffusion coefficient. It is the CPU/GPU-portable analogue of Open3D's
    ``TriangleMesh.filter_smooth_simple``.

    Isolated vertices (degree 0) are left unchanged.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    iterations
        Number of smoothing passes.
    laplacian_operator
        Optional precomputed row-stochastic operator (see
        [`laplacian`][triwarp.laplacian.laplacian]). Autogenerated (uniform weights) when ``None``.

    Returns
    -------
    wp.array[wp.vec3]
        Smoothed ``(n_vertices,)`` vertex positions on ``vertices.device``.

    Raises
    ------
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian]
    [`filter_taubin`][triwarp.smoothing.filter_taubin]
    """
    require_same_device(vertices=vertices, faces=faces)
    device = vertices.device
    n = int(vertices.shape[0])
    if n == 0 or iterations == 0:
        return wp.clone(vertices)

    # Symmetric adjacency (undirected 1-ring) so boundary vertices average over all their
    # neighbors, matching Open3D's ``adjacency_list``; the directed default is asymmetric there.
    operator = _resolved_operator(vertices, faces, laplacian_operator, symmetric=True)
    positions = _as_vec3d(vertices)
    nxt = wp.empty(n, dtype=wp.vec3d, device=device)
    for _ in range(iterations):
        # The kernel reads its own degree off the CSR row it is already walking:
        # degree(i) = offsets[i + 1] - offsets[i].
        wp.launch(
            kernel_smoothing.neighborhood_average_pass,
            dim=n,
            inputs=[operator.offsets, operator.columns, operator.values, positions],
            outputs=[nxt],
            device=device,
        )
        positions, nxt = nxt, positions

    return _as_vec3(positions)

filter_mut_dif_laplacian 

filter_mut_dif_laplacian(
    vertices: array[vec3],
    faces: array[int32],
    lamb: float = 0.5,
    iterations: int = 10,
    volume_constraint: bool = True,
    laplacian_operator: BsrMatrix[float32] | None = None,
    *,
    face_normals: array[vec3] | None = None,
    face_areas: array[float32] | None = None,
) -> wp.array[wp.vec3]

Laplacian smoothing with a mutable diffusion coefficient (Barroqueiro et al.).

Each pass runs a Laplacian diffusion step whose speed is adapted per vertex: the coefficient lamber is scaled by how much of the vertex's Laplacian residual lies along its normal, so flat regions diffuse quickly while feature regions diffuse slowly. When volume_constraint is set, the mesh is inflated along its vertex normals (rather than isotropically rescaled) to counteract shrinkage. Mirrors trimesh.smoothing.filter_mut_dif_laplacian but returns a new vertex array instead of mutating the mesh in place.

Following the reference, the vertex normals (and the finite-difference step eps) are computed once from the input mesh and held constant across all passes.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
lamb float

Diffusion speed. 0 leaves the mesh unchanged; larger values smooth faster.

 0.5
iterations int

Number of smoothing passes.

 10
volume_constraint bool

If True inflate the mesh along its vertex normals after each pass to preserve its initial volume (meaningful only for watertight meshes).

 True
laplacian_operator BsrMatrix[float32] | None

Optional precomputed row-stochastic operator (see laplacian). Autogenerated (uniform weights) when None.

 None
face_normals array[vec3] | None

Optional length-n_faces unit face normals and matching areas from face_normals_and_areas; recomputed together when either is None. Trimesh.face_normals and Trimesh.face_areas cache the pair.

 None
face_areas array[float32] | None

See face_normals.

 None

Returns:

Type Description
array[vec3]

Smoothed (n_vertices,) vertex positions on vertices.device.

Raises:

Type Description
RuntimeError

If vertices, faces, face_normals and face_areas are not all on one device.
See Also

filter_laplacian filter_humphrey filter_taubin trimesh.smoothing.filter_mut_dif_laplacian

Source code in triwarp/smoothing.py 
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
def filter_mut_dif_laplacian(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    lamb: float = 0.5,
    iterations: int = 10,
    volume_constraint: bool = True,
    laplacian_operator: wps.BsrMatrix[wp.float32] | None = None,
    *,
    face_normals: wp.array[wp.vec3] | None = None,
    face_areas: wp.array[wp.float32] | None = None,
) -> wp.array[wp.vec3]:
    """
    Laplacian smoothing with a mutable diffusion coefficient (Barroqueiro et al.).

    Each pass runs a Laplacian diffusion step whose speed is adapted **per vertex**: the
    coefficient ``lamber`` is scaled by how much of the vertex's Laplacian residual lies along
    its normal, so flat regions diffuse quickly while feature regions diffuse slowly. When
    ``volume_constraint`` is set, the mesh is inflated along its vertex normals (rather than
    isotropically rescaled) to counteract shrinkage. Mirrors
    [`trimesh.smoothing.filter_mut_dif_laplacian`][] but returns a new vertex array instead of
    mutating the mesh in place.

    Following the reference, the vertex normals (and the finite-difference step ``eps``) are
    computed once from the input mesh and held constant across all passes.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    lamb
        Diffusion speed. ``0`` leaves the mesh unchanged; larger values smooth faster.
    iterations
        Number of smoothing passes.
    volume_constraint
        If ``True`` inflate the mesh along its vertex normals after each pass to preserve its
        initial volume (meaningful only for watertight meshes).
    laplacian_operator
        Optional precomputed row-stochastic operator (see
        [`laplacian`][triwarp.laplacian.laplacian]). Autogenerated (uniform weights) when ``None``.
    face_normals
        Optional length-``n_faces`` unit face normals and matching areas from
        [`face_normals_and_areas`][triwarp.triangles.face_normals_and_areas]; recomputed together
        when either is ``None``. [`Trimesh.face_normals`][triwarp.mesh.Trimesh.face_normals] and
        [`Trimesh.face_areas`][triwarp.mesh.Trimesh.face_areas] cache the pair.
    face_areas
        See ``face_normals``.

    Returns
    -------
    wp.array[wp.vec3]
        Smoothed ``(n_vertices,)`` vertex positions on ``vertices.device``.

    Raises
    ------
    RuntimeError
        If ``vertices``, ``faces``, ``face_normals`` and ``face_areas`` are not all on one device.

    See Also
    --------
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian]
    [`filter_humphrey`][triwarp.smoothing.filter_humphrey]
    [`filter_taubin`][triwarp.smoothing.filter_taubin]
    [`trimesh.smoothing.filter_mut_dif_laplacian`][]
    """
    require_same_device(
        vertices=vertices, faces=faces, face_normals=face_normals, face_areas=face_areas
    )
    device = vertices.device
    n = int(vertices.shape[0])
    if n == 0 or iterations == 0:
        return wp.clone(vertices)

    operator = _resolved_operator(vertices, faces, laplacian_operator)
    positions = _as_vec3d(vertices)

    # Vertex normals and eps are computed once from the input mesh and reused every pass, matching
    # the trimesh reference (which reads normals off the un-mutated mesh inside its loop).
    if face_normals is None or face_areas is None:
        face_normals, face_areas = face_normals_and_areas(vertices, faces)
    normals = mean_vertex_normals(n, faces, face_normals)
    vol_ini = tw.measures.volume(positions, faces) if volume_constraint else 0.0
    eps = 0.01 * float(tw.reduce.max(face_areas)) ** 0.5 if volume_constraint else 0.0

    lv = wp.empty(n, dtype=wp.vec3d, device=device)
    adil = wp.empty(n, dtype=wp.float64, device=device)
    adil_sum = wp.zeros(1, dtype=wp.float64, device=device)
    nxt = wp.empty(n, dtype=wp.vec3d, device=device)
    probe = wp.empty(n, dtype=wp.vec3d, device=device) if volume_constraint else None
    slope = 0.0
    inv_n = wp.float64(1.0 / n)
    n_blocks = kernel_reduce.blocks_1d(n)
    for index in range(iterations):
        # The mean diffusion coefficient is reduced on device and consumed by the step kernel
        # directly, so the loop body issues no host synchronisation.
        wp.launch(
            kernel_smoothing.mut_dif_adil_pass,
            dim=n,
            inputs=[operator.offsets, operator.columns, operator.values, positions, normals],
            outputs=[lv, adil],
            device=device,
        )
        adil_sum.zero_()
        wp.launch_tiled(
            kernel_reduce.SUM1D_TILED[wp.float64],
            dim=[n_blocks],
            inputs=[adil, adil_sum],
            block_dim=TILE_1D,
            device=device,
        )
        wp.launch(
            kernel_smoothing.mut_dif_step_scaled,
            dim=n,
            inputs=[positions, lv, adil, adil_sum, inv_n, wp.float64(lamb)],
            outputs=[nxt],
            device=device,
        )
        positions, nxt = nxt, positions
        if volume_constraint:
            vol = tw.measures.volume(positions, faces)
            if index == 0:
                wp.map(
                    kernel_smoothing.add_scaled_normal,
                    positions,
                    normals,
                    wp.float64(eps),
                    out=probe,
                )
                vol2 = tw.measures.volume(probe, faces)
                slope = eps / (vol2 - vol) if vol2 != vol else 0.0
            wp.map(
                kernel_smoothing.add_scaled_normal,
                positions,
                normals,
                wp.float64(slope * (vol_ini - vol)),
                out=positions,
            )

    return _as_vec3(positions)

filter_implicit_fairing 

filter_implicit_fairing(
    vertices: array[vec3],
    faces: array[int32],
    lamb: float = 0.1,
    iterations: int = 10,
    *,
    pin_boundary: bool = True,
) -> wp.array[wp.vec3]

Implicit fairing with the cotangent Laplace-Beltrami operator (Desbrun et al.).

Solves the backward-Euler diffusion (M - lamb L) V' = M V each pass, where L is the cotangent stiffness matrix (cotmatrix) and M the barycentric lumped mass matrix (mass_matrix), both recomputed from the current geometry. This is the most accurate of the smoothing filters (curvature flow; libigl tutorial 205). The solve uses float64 throughout.

On a mesh with an open boundary the flow needs a boundary condition, which is what pin_boundary supplies: without one the unconstrained boundary is pulled inward, collapsing the triangles there, and since the system is solved by conjugate gradient (the only Warp solver available) further passes diverge rather than degrade gracefully. Pinning turns the pass into a Dirichlet problem over the interior, which is the formulation Desbrun et al. give for a surface with boundary and is well posed for any number of passes.

A closed mesh has no boundary, so pin_boundary costs nothing and changes nothing there -- the solve is over every vertex either way, and curvature flow on a closed surface still shrinks it (a sphere contracts toward a point, which is the method working, not failing).

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
lamb float

Diffusion time step. Larger values smooth faster; the implicit solve is unconditionally stable.

 0.1
iterations int

Number of smoothing passes.

 10
pin_boundary bool

Hold boundary vertices at their input positions and solve only for the interior. A no-op on a closed mesh. Set False for the unconstrained flow, which is free to shrink the boundary — expect solve_spd to warn once it stops converging.

 True

Returns:

Type Description
array[vec3]

Smoothed (n_vertices,) vertex positions on vertices.device.

Raises:

Type Description
RuntimeError

If vertices and faces are not all on one device.
See Also

filter_laplacian cotmatrix mass_matrix

Source code in triwarp/smoothing.py 
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
def filter_implicit_fairing(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    lamb: float = 0.1,
    iterations: int = 10,
    *,
    pin_boundary: bool = True,
) -> wp.array[wp.vec3]:
    """
    Implicit fairing with the cotangent Laplace-Beltrami operator (Desbrun et al.).

    Solves the backward-Euler diffusion ``(M - lamb L) V' = M V`` each pass, where ``L`` is the
    cotangent stiffness matrix ([`cotmatrix`][triwarp.laplacian.cotmatrix]) and ``M`` the
    barycentric lumped mass matrix ([`mass_matrix`][triwarp.laplacian.mass_matrix]), both recomputed
    from the current geometry. This is the most accurate of the smoothing filters (curvature flow;
    libigl tutorial 205). The solve uses ``float64`` throughout.

    On a mesh with an open boundary the flow needs a boundary condition, which is what
    ``pin_boundary`` supplies: without one the unconstrained boundary is pulled inward, collapsing
    the triangles there, and since the system is solved by conjugate gradient (the only Warp solver
    available) further passes diverge rather than degrade gracefully. Pinning turns the pass into a
    Dirichlet problem over the interior, which is the formulation Desbrun et al. give for a surface
    with boundary and is well posed for any number of passes.

    A closed mesh has no boundary, so ``pin_boundary`` costs nothing and changes nothing there --
    the solve is over every vertex either way, and curvature flow on a closed surface still shrinks
    it (a sphere contracts toward a point, which is the method working, not failing).

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    lamb
        Diffusion time step. Larger values smooth faster; the implicit solve is unconditionally
        stable.
    iterations
        Number of smoothing passes.
    pin_boundary
        Hold boundary vertices at their input positions and solve only for the interior. A no-op on
        a closed mesh. Set ``False`` for the unconstrained flow, which is free to shrink the
        boundary — expect [`solve_spd`][triwarp.linalg.solve_spd] to warn once it stops converging.

    Returns
    -------
    wp.array[wp.vec3]
        Smoothed ``(n_vertices,)`` vertex positions on ``vertices.device``.

    Raises
    ------
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian]
    [`cotmatrix`][triwarp.laplacian.cotmatrix]
    [`mass_matrix`][triwarp.laplacian.mass_matrix]
    """
    require_same_device(vertices=vertices, faces=faces)
    device = vertices.device
    n = int(vertices.shape[0])
    n_faces = int(faces.shape[0]) // 3
    if n == 0 or n_faces == 0 or iterations == 0:
        return wp.clone(vertices)

    positions = _as_vec3d(vertices)
    components = _component_columns(n, device)
    rhs = _component_columns(n, device)
    solutions = _component_columns(n, device)
    component_rows = [components[column] for column in range(3)]

    # Boundary topology is fixed for the whole flow, so the partition is built once. ``None`` means
    # "solve over every vertex" -- either the caller asked for the unconstrained flow, or the mesh
    # is closed and there is nothing to pin.
    dirichlet = _dirichlet_state(vertices, faces, positions, n, device) if pin_boundary else None
    if dirichlet is not None and dirichlet.n_free == 0:
        # Every vertex is pinned, so the flow has no unknown to move and the pass is the identity.
        return _as_vec3(positions)

    for _ in range(iterations):
        current = _as_vec3(positions)
        # Rebuilt every iteration on purpose: the cotangent weights depend on ``current``, which
        # the fairing step moves, so implicit fairing must re-linearise on the moving surface.
        cot_entries = laplacian.cotmatrix_entries(current, faces)
        stiffness = laplacian.cotmatrix(current, faces, cot_entries=cot_entries, dtype=wp.float64)

        mass = laplacian.mass_matrix_entries(current, faces, dtype=wp.float64)

        # Right-hand side b = M V, formed before ``bsr_axpy`` mutates the mass matrix.
        wp.map(kernel_smoothing.extract_components, positions, out=component_rows)
        for column in range(3):
            wp.map(wp.mul, mass, components[column], out=rhs[column])

        # A = M - lamb L (SPD: L has a negative diagonal, so subtracting it adds to the diagonal).
        system = wps.bsr_axpy(x=stiffness, y=wps.bsr_diag(diag=mass), alpha=-float(lamb), beta=1.0)

        if dirichlet is None:
            # Seed CG with the current positions, not with ``b = M V``: an unreferenced vertex has
            # an all-zero row and a zero right-hand side, so CG never writes its entry and it would
            # keep whatever the seed left there. One batched solve advances all three columns
            # together; the operator is rebuilt every pass, so there is no state to hoist.
            wp.copy(solutions, components)
            twl.solve_spd_columns(system, rhs, solutions, tol=twl.CG_TOLERANCE, maxiter=10 * n)
            wp.map(
                kernel_smoothing.combine_components,
                solutions[0],
                solutions[1],
                solutions[2],
                out=positions,
            )
            continue

        # Dirichlet pass: eliminate the pinned rows and columns, then solve over the interior.
        # ``assemble_interior_system`` supplies ``-A_ub x_b``; the linear term ``(M V)_u`` is added
        # on top, since that helper eliminates a quadratic form which has none of its own.
        interior_system, interior_rhs = twl.assemble_interior_system(
            system, dirichlet.fixed_mask, dirichlet.free_map, dirichlet.pinned, dirichlet.n_free
        )
        wp.launch(
            kernel_smoothing.add_interior_mass_rhs,
            dim=n,
            inputs=[dirichlet.fixed_mask, dirichlet.free_map, mass, positions, interior_rhs],
            device=device,
        )
        solution_2d = dirichlet.solution
        # Seed CG with the free vertices' current positions, exactly as the unreduced branch above
        # does -- the right-hand side is not a position, and an unreferenced vertex whose row and
        # right-hand side are both zero would be left wherever the seed put it.
        wp.launch(
            kernel_smoothing.gather_free_positions_2d,
            dim=n,
            inputs=[dirichlet.fixed_mask, dirichlet.free_map, positions, solution_2d],
            device=device,
        )
        twl.solve_spd_columns(
            interior_system,
            interior_rhs,
            solution_2d,
            tol=twl.CG_TOLERANCE,
            maxiter=10 * dirichlet.n_free,
        )
        wp.launch(
            kernel_smoothing.scatter_free_positions,
            dim=n,
            inputs=[
                dirichlet.fixed_mask,
                dirichlet.free_map,
                solution_2d[0],
                solution_2d[1],
                solution_2d[2],
                positions,
            ],
            device=device,
        )

    return _as_vec3(positions)

smooth_region_fixed_rim 

smooth_region_fixed_rim(
    vertices: array[vec3],
    faces: array[int32],
    free_mask: array[bool],
    stabilizer: float = 0.0,
) -> wp.array[wp.vec3]

Reposition a free vertex region as the umbrella-Laplacian solution with a sharp fixed boundary.

The free vertices (free_mask) are moved to the solution of the graph-Laplacian Dirichlet system (D - W) x = b (unit edge weights), where the fixed one-ring neighbours are folded into b, so the region rim stays sharp (a hard C⁰ constraint). Fixed vertices keep their positions.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) vertex positions.

 required
faces array[int32]

Length-3 * n_faces wp.int32 flat triangle index buffer.

 required
free_mask array[bool]

Length-n_vertices wp.bool mask of the vertices to reposition.

 required
stabilizer float

Optional attraction to the original position (adds to the diagonal), which also makes the system solvable for a fully-free component. Defaults to 0.

 0.0

Returns:

Type Description
array[vec3]

New vertex positions on vertices.device (a copy; fixed vertices unchanged).

Raises:

Type Description
RuntimeError

If vertices, faces and free_mask are not all on one device.
See Also

smooth_region smooth_region_boundary The third case: smooths the rim's own path rather than the surface on either side of it. fill_smooth

Notes

The system is SPD only when no free connected component is entirely free (or stabilizer > 0); the hole-filling pipeline guarantees this because the patch rim is always fixed.

A free vertex that no face refers to contributes no row, so the solve does not constrain it and it is returned where it was. That is what the solution seed decides rather than the system: it is seeded with the current positions, not with zeros.

Source code in triwarp/smoothing.py 
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
def smooth_region_fixed_rim(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    free_mask: wp.array[wp.bool],
    stabilizer: float = 0.0,
) -> wp.array[wp.vec3]:
    """
    Reposition a free vertex region as the umbrella-Laplacian solution with a sharp fixed boundary.

    The free vertices (``free_mask``) are moved to the solution of the graph-Laplacian Dirichlet
    system ``(D - W) x = b`` (unit edge weights), where the fixed one-ring neighbours are folded
    into ``b``, so the region rim stays sharp
    (a hard C⁰ constraint). Fixed vertices keep their positions.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer.
    free_mask
        Length-``n_vertices`` ``wp.bool`` mask of the vertices to reposition.
    stabilizer
        Optional attraction to the original position (adds to the diagonal), which also makes the
        system solvable for a fully-free component. Defaults to ``0``.

    Returns
    -------
    wp.array[wp.vec3]
        New vertex positions on ``vertices.device`` (a copy; fixed vertices unchanged).

    Raises
    ------
    RuntimeError
        If ``vertices``, ``faces`` and ``free_mask`` are not all on one device.

    See Also
    --------
    [`smooth_region`][triwarp.smoothing.smooth_region]
    [`smooth_region_boundary`][triwarp.smoothing.smooth_region_boundary]
        The third case: smooths the rim's own path rather than the surface on either side of it.
    [`fill_smooth`][triwarp.holes.fill_smooth]

    Notes
    -----
    The system is SPD only when no free connected component is entirely free (or ``stabilizer >
    0``); the hole-filling pipeline guarantees this because the patch rim is always fixed.

    A free vertex that no face refers to contributes no row, so the solve does not constrain it
    and it is returned where it was. That is what the ``solution`` seed decides rather than the
    system: it is seeded with the current positions, not with zeros.
    """
    require_same_device(vertices=vertices, faces=faces, free_mask=free_mask)
    device = vertices.device
    n = int(vertices.shape[0])
    out = wp.clone(vertices)
    if int(faces.shape[0]) == 0 or n == 0:
        return out
    free_map, n_free = tw.array.mask_to_compact_ranks(free_mask)
    if n_free == 0:
        return out

    weight_matrix = _edge_weight_matrix(vertices, faces, "unit")
    nnz = int(weight_matrix.nnz)
    size = nnz + n
    # ``dirichlet_system_triplets`` emits conditionally, so most slots stay unwritten. Padding must
    # be a *hole*, not a value: a ``(0, 0, 0.0)`` triplet is a harmless structural zero, but every
    # one of them accumulates onto entry ``(0, 0)`` and ``bsr_from_triplets``' accumulation atomic
    # serializes them. An out-of-range index is dropped instead -- ``n_free`` is one past the last
    # row and column of the ``(n_free, n_free)`` system.
    out_rows = wp.full(size, n_free, dtype=wp.int32, device=device)
    out_cols = wp.full(size, n_free, dtype=wp.int32, device=device)
    out_vals = wp.zeros(size, dtype=wp.float64, device=device)
    # One contiguous (3, n_free) right-hand side: its rows are contiguous 1-D views, so the
    # assembly kernel writes them directly and the three columns solve in one batched CG.
    rhs = wp.zeros((3, n_free), dtype=wp.float64, device=device)
    wp.launch(
        kernel_smoothing.dirichlet_system_triplets,
        dim=n,
        inputs=[
            weight_matrix.offsets,
            weight_matrix.columns,
            weight_matrix.values,
            free_mask,
            free_map,
            vertices,
            wp.float64(stabilizer),
            out_rows,
            out_cols,
            out_vals,
            rhs[0],
            rhs[1],
            rhs[2],
        ],
        device=device,
    )
    system = wps.bsr_from_triplets(
        n_free, n_free, out_rows, out_cols, out_vals, prune_numerical_zeros=False
    )
    sol = _free_positions(vertices, free_mask, free_map, n_free)
    twl.solve_spd_columns(
        system,
        twt.as_array2d(rhs, wp.float64),
        twt.as_array2d(sol, wp.float64),
        tol=twl.CG_TOLERANCE,
        maxiter=10 * n_free,
    )
    wp.launch(
        kernel_smoothing.scatter_free_solution,
        dim=n,
        inputs=[free_mask, free_map, sol[0], sol[1], sol[2], out],
        device=device,
    )
    return out

smooth_region 

smooth_region(
    vertices: array[vec3],
    faces: array[int32],
    free_mask: array[bool],
    edge_weights: str = "cotan",
) -> wp.array[wp.vec3]

Reposition a free vertex region so the surface is smooth across the region boundary too.

Every vertex in the region and its first fixed ring contributes the umbrella equation p_v = Σ_d (w_vd / ΣW) p_d; free vertices are unknowns and fixed neighbours move to the right-hand side. The normal equations (MᵀM) x = Mᵀ b are solved per coordinate, giving a patch that is smooth (C¹) across the region rim, unlike the sharp-rim smooth_region_fixed_rim.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) vertex positions.

 required
faces array[int32]

Length-3 * n_faces wp.int32 flat triangle index buffer.

 required
free_mask array[bool]

Length-n_vertices wp.bool mask of the vertices to reposition.

 required
edge_weights str

"cotan" (default, clamped cotangent weights) or "unit" (uniform).

 'cotan'

Returns:

Type Description
array[vec3]

New vertex positions on vertices.device (a copy; fixed vertices unchanged).

Raises:

Type Description
ValueError

If edge_weights is not "cotan" or "unit".
RuntimeError

If vertices, faces and free_mask are not all on one device.
See Also

smooth_region_fixed_rim smooth_region_boundary The third case: smooths the rim's own path rather than the surface on either side of it. fill_smooth

Notes

A free vertex that no face refers to contributes no row, so the least-squares system does not constrain it and it is returned where it was. That is what the solution seed decides rather than the system: it is seeded with the current positions, not with zeros.

Source code in triwarp/smoothing.py 
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
def smooth_region(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    free_mask: wp.array[wp.bool],
    edge_weights: str = "cotan",
) -> wp.array[wp.vec3]:
    """
    Reposition a free vertex region so the surface is smooth across the region boundary too.

    Every vertex in the region and its first fixed ring contributes the umbrella equation
    ``p_v = Σ_d (w_vd / ΣW) p_d``; free vertices are unknowns and fixed neighbours move to the
    right-hand side. The normal equations ``(MᵀM) x = Mᵀ b`` are solved per
    coordinate, giving a patch that is smooth (C¹) *across* the region rim, unlike the sharp-rim
    [`smooth_region_fixed_rim`][triwarp.smoothing.smooth_region_fixed_rim].

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer.
    free_mask
        Length-``n_vertices`` ``wp.bool`` mask of the vertices to reposition.
    edge_weights
        ``"cotan"`` (default, clamped cotangent weights) or ``"unit"`` (uniform).

    Returns
    -------
    wp.array[wp.vec3]
        New vertex positions on ``vertices.device`` (a copy; fixed vertices unchanged).

    Raises
    ------
    ValueError
        If ``edge_weights`` is not ``"cotan"`` or ``"unit"``.
    RuntimeError
        If ``vertices``, ``faces`` and ``free_mask`` are not all on one device.

    See Also
    --------
    [`smooth_region_fixed_rim`][triwarp.smoothing.smooth_region_fixed_rim]
    [`smooth_region_boundary`][triwarp.smoothing.smooth_region_boundary]
        The third case: smooths the rim's own path rather than the surface on either side of it.
    [`fill_smooth`][triwarp.holes.fill_smooth]

    Notes
    -----
    A free vertex that no face refers to contributes no row, so the least-squares system does not
    constrain it and it is returned where it was. That is what the ``solution`` seed decides rather
    than the system: it is seeded with the current positions, not with zeros.
    """
    require_same_device(vertices=vertices, faces=faces, free_mask=free_mask)
    device = vertices.device
    n = int(vertices.shape[0])
    out = wp.clone(vertices)
    if int(faces.shape[0]) == 0 or n == 0:
        return out
    free_map, n_free = tw.array.mask_to_compact_ranks(free_mask)
    if n_free == 0:
        return out

    row_mask = tw.selection.expand_vertex_mask(faces, free_mask, 1)
    row_map, n_rows = tw.array.mask_to_compact_ranks(row_mask)
    weight_matrix = _edge_weight_matrix(vertices, faces, edge_weights)
    nnz = int(weight_matrix.nnz)
    size = nnz + n
    # ``laplacian_ls_triplets`` emits conditionally, so most slots stay unwritten. Padding must be a
    # *hole*, not a value: a ``(0, 0, 0.0)`` triplet is a harmless structural zero, but every one of
    # them accumulates onto entry ``(0, 0)`` and ``bsr_from_triplets``' accumulation atomic
    # serializes them. Both arrays are filled one past their own extent, so the padding is out of
    # range as the *row* index of ``M`` (``rows``) and of ``M^T`` (``cols``) alike.
    rows = wp.full(size, n_rows, dtype=wp.int32, device=device)
    cols = wp.full(size, n_free, dtype=wp.int32, device=device)
    vals = wp.zeros(size, dtype=wp.float64, device=device)
    rhs_x = wp.zeros(n_rows, dtype=wp.float64, device=device)
    rhs_y = wp.zeros(n_rows, dtype=wp.float64, device=device)
    rhs_z = wp.zeros(n_rows, dtype=wp.float64, device=device)
    wp.launch(
        kernel_smoothing.laplacian_ls_triplets,
        dim=n,
        inputs=[
            weight_matrix.offsets,
            weight_matrix.columns,
            weight_matrix.values,
            free_mask,
            row_mask,
            free_map,
            row_map,
            vertices,
            rows,
            cols,
            vals,
            rhs_x,
            rhs_y,
            rhs_z,
        ],
        device=device,
    )
    # M is (n_rows x n_free) and A = M^T M is SPD.
    #
    # ``m_matrix``'s ``nnz`` is a padded *capacity*, not an exact count: ``bsr_from_triplets``
    # leaves the field at the triplet count it was handed, which is the ``size`` above --
    # overwhelmingly padding here, since the emit is conditional and the tail is deliberately out
    # of range. That is harmless to ``bsr_mm`` and ``bsr_transposed``, which treat the field as a
    # bound rather than as a count, but nothing downstream may be sized off it (section 3.7).
    #
    # The transpose used to be a second ``bsr_from_triplets`` over the same triplets with the two
    # index arrays swapped, plus three ``wp.clone``s to feed it. ``bsr_transposed`` is exact on an
    # operand of this shape and says in one call what the swap said in five lines.
    #
    # **Not a speed change, and the number is here so it is not re-proposed as one.** Interleaved
    # against the swapped-triplet build at the benchmark's own operating point (the top quarter by
    # z of an icosphere, min of 4 rounds of 10), it is 1.16x / 1.14x / 1.09x on the step at 2 562 /
    # 10 242 / 40 962 vertices -- which is 0.41% / 0.19% / 0.05% of the whole call, a share that
    # *falls* as the mesh grows and is therefore a decline by section 9's rule. What it buys is one
    # concept fewer and three fewer buffers; the assembled system is bit-identical at all three
    # sizes (max |delta| exactly 0.0 over the sorted values, same nnz).
    m_matrix = wps.bsr_from_triplets(n_rows, n_free, rows, cols, vals, prune_numerical_zeros=False)
    mt_matrix = wps.bsr_transposed(m_matrix)
    system = wps.bsr_mm(mt_matrix, m_matrix)
    # A^T b straight into the rows of one contiguous buffer, so the three columns batch.
    atb = wp.zeros((3, n_free), dtype=wp.float64, device=device)
    for column, component in enumerate((rhs_x, rhs_y, rhs_z)):
        wps.bsr_mv(mt_matrix, component, atb[column], alpha=1.0, beta=0.0)
    sol = _free_positions(vertices, free_mask, free_map, n_free)
    # The normal equations here are the worst-conditioned system this package solves, so a
    # multigrid hierarchy usually pays for its own setup. ``"auto"`` rather than an unconditional
    # ``"multigrid"`` because this same function is also called on small hole patches by
    # ``refine_and_smooth_region``, where the hierarchy's setup cost is not repaid. What separates a
    # patch from a region is not size or iteration count alone but the operator's off-diagonal
    # dominance, which is what ``"auto"`` gates on (``linalg.CG_MULTIGRID_DOMINANCE``).
    twl.solve_spd_columns(
        system,
        twt.as_array2d(atb, wp.float64),
        twt.as_array2d(sol, wp.float64),
        tol=twl.CG_TOLERANCE,
        maxiter=10 * n_free,
        preconditioner="auto",
    )
    wp.launch(
        kernel_smoothing.scatter_free_solution,
        dim=n,
        inputs=[free_mask, free_map, sol[0], sol[1], sol[2], out],
        device=device,
    )
    return out

refine_and_smooth_region 

refine_and_smooth_region(
    vertices: array[vec3],
    faces: array[int32],
    n_vertices_before: int,
    patch_face_mask: array[bool],
    max_edge: float,
    max_edge_splits: int,
    max_angle_change_after_flip: float,
    smooth_curvature: bool,
    smooth_boundary: bool,
    natural_smooth: bool,
    edge_weights: str,
    *,
    refine: Literal["max_edge", "density"] = "max_edge",
) -> tuple[
    wp.array[wp.vec3], wp.array[wp.int32], wp.array[wp.bool]
]

Subdivide the patch and smooth its new vertices.

Shared finisher of fill_smooth and stitch_smooth.

refine selects the subdivision criterion, and the two are genuinely different questions rather than two tunings of one: "max_edge" bisects region edges longer than max_edge (remesh.subdivide_region_to_size), and "density" splits region triangles at their centroid while their sampling is coarser than the surrounding mesh's -- Liepa's criterion, in remesh.refine_region_to_density. max_edge and max_edge_splits are ignored under "density", which takes its target from the mesh instead of from an argument.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions, patch included.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer, patch included.

 required
n_vertices_before int

Vertex count before the patch was appended. Everything at or above this index is a patch vertex, which is how the smooth step finds the vertices it may move.

 required
patch_face_mask array[bool]

(n_faces,) flag marking the patch faces to refine and smooth.

 required
max_edge float

Target edge length for refine="max_edge". Ignored under "density".

 required
max_edge_splits int

Soft cap on the number of edge splits during subdivision. Ignored under "density".

 required
max_angle_change_after_flip float

Dihedral-angle-change gate for the Delaunay flip pass, in radians.

 required
smooth_curvature bool

When True, smooth the new patch vertices after subdivision. When False the function returns straight after refining.

 required
smooth_boundary bool

When True, let the patch rim move as well as its interior.

 required
natural_smooth bool

When True, additionally grow a collar around the patch and smooth it so the patch blends into the surrounding surface.

 required
edge_weights str

Laplacian edge weights for the smooth solve: "cotan" or "unit".

 required
refine Literal['max_edge', 'density']

Subdivision criterion, as above.

 'max_edge'

Returns:

Name Type Description
vertices array[vec3]

(n_vertices,) positions, the input's plus whatever subdivision appended, with the patch smoothed when smooth_curvature.
faces array[int32]

Flat 3 * n_faces triangle index buffer of the refined mesh.
patch_face_mask array[bool]

The patch mask carried through the refinement, one flag per face of the returned buffer -- longer than the input mask, since subdivision adds faces.

Raises:

Type Description
ValueError

If refine is neither "max_edge" nor "density".
RuntimeError

If vertices, faces and patch_face_mask are not all on one device.
Source code in triwarp/smoothing.py 
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
1934
1935
1936
1937
1938
1939
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
def refine_and_smooth_region(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    n_vertices_before: int,
    patch_face_mask: wp.array[wp.bool],
    max_edge: float,
    max_edge_splits: int,
    max_angle_change_after_flip: float,
    smooth_curvature: bool,
    smooth_boundary: bool,
    natural_smooth: bool,
    edge_weights: str,
    *,
    refine: Literal["max_edge", "density"] = "max_edge",
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32], wp.array[wp.bool]]:
    """
    Subdivide the patch and smooth its new vertices.

    Shared finisher of [`fill_smooth`][triwarp.holes.fill_smooth] and
    [`stitch_smooth`][triwarp.holes.stitch_smooth].

    ``refine`` selects the subdivision criterion, and the two are genuinely different questions
    rather than two tunings of one: ``"max_edge"`` bisects region edges longer than ``max_edge``
    ([`remesh.subdivide_region_to_size`][triwarp.remesh.subdivide_region_to_size]), and
    ``"density"`` splits region triangles at their centroid while their sampling is coarser than
    the surrounding mesh's -- Liepa's criterion, in
    [`remesh.refine_region_to_density`][triwarp.remesh.refine_region_to_density]. ``max_edge`` and
    ``max_edge_splits`` are ignored under ``"density"``, which takes its target from the mesh
    instead of from an argument.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions, patch included.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer, patch included.
    n_vertices_before
        Vertex count before the patch was appended. Everything at or above this index is a patch
        vertex, which is how the smooth step finds the vertices it may move.
    patch_face_mask
        ``(n_faces,)`` flag marking the patch faces to refine and smooth.
    max_edge
        Target edge length for ``refine="max_edge"``. Ignored under ``"density"``.
    max_edge_splits
        Soft cap on the number of edge splits during subdivision. Ignored under ``"density"``.
    max_angle_change_after_flip
        Dihedral-angle-change gate for the Delaunay flip pass, in radians.
    smooth_curvature
        When ``True``, smooth the new patch vertices after subdivision. When ``False`` the function
        returns straight after refining.
    smooth_boundary
        When ``True``, let the patch rim move as well as its interior.
    natural_smooth
        When ``True``, additionally grow a collar around the patch and smooth it so the patch blends
        into the surrounding surface.
    edge_weights
        Laplacian edge weights for the smooth solve: ``"cotan"`` or ``"unit"``.
    refine
        Subdivision criterion, as above.

    Returns
    -------
    vertices : wp.array[wp.vec3]
        ``(n_vertices,)`` positions, the input's plus whatever subdivision appended, with the patch
        smoothed when ``smooth_curvature``.
    faces : wp.array[wp.int32]
        Flat ``3 * n_faces`` triangle index buffer of the refined mesh.
    patch_face_mask : wp.array[wp.bool]
        The patch mask carried through the refinement, one flag per face of the *returned* buffer --
        longer than the input mask, since subdivision adds faces.

    Raises
    ------
    ValueError
        If ``refine`` is neither ``"max_edge"`` nor ``"density"``.
    RuntimeError
        If ``vertices``, ``faces`` and ``patch_face_mask`` are not all on one device.
    """
    require_same_device(vertices=vertices, faces=faces, patch_face_mask=patch_face_mask)
    device = faces.device
    if refine == "density":
        vertices, faces, patch_face_mask = tw.remesh.refine_region_to_density(
            vertices, faces, patch_face_mask, max_angle_change=max_angle_change_after_flip
        )
    elif refine == "max_edge":
        vertices, faces, patch_face_mask = tw.remesh.subdivide_region_to_size(
            vertices,
            faces,
            patch_face_mask,
            max_edge=max_edge,
            max_splits=max_edge_splits,
            max_angle_change=max_angle_change_after_flip,
        )
    else:
        raise ValueError(f"unknown refine {refine!r}, expected 'max_edge' or 'density'")
    if not smooth_curvature:
        return vertices, faces, patch_face_mask

    n = int(vertices.shape[0])
    # New (interior patch) vertices are the tail appended by subdivision, minus mesh-boundary verts.
    new_verts = wp.zeros(n, dtype=wp.bool, device=device)
    # Warp rejects a zero-length slice, and subdivision may have added no vertices at all.
    if n > n_vertices_before:
        new_verts[n_vertices_before:].fill_(True)
    bd_mask = _boundary_verts_mask(vertices, faces)
    free = wp.empty(n, dtype=wp.bool, device=device)
    wp.map(kernel_array.mask_and_not, new_verts, bd_mask, out=free)

    vertices = smooth_region_fixed_rim(vertices, faces, free)
    if smooth_boundary:
        vertices = smooth_region(vertices, faces, free, edge_weights)

    if natural_smooth:
        edges_bd = tw.selection.region_boundary_edges(faces, patch_face_mask, n_vertices=n)
        endpoints = wp.clone(twt.as_array2d(edges_bd, wp.int32).reshape(-1))
        incident = tw.array.indices_to_mask(endpoints, n, device=device)
        incident = tw.selection.expand_vertex_mask(faces, incident, 5)
        incident = tw.selection.shrink_vertex_mask(faces, incident, 2)
        incident = tw.selection.exclude_fully_selected_components(faces, incident, n)
        # A one-byte-per-vertex device reduction is cheap enough here to prefer over a host
        # readback.
        if tw.reduce.any(incident):
            bd_mask = _boundary_verts_mask(vertices, faces)
            free2 = wp.empty(n, dtype=wp.bool, device=device)
            wp.map(kernel_array.mask_and_not, incident, bd_mask, out=free2)
            vertices = smooth_region_fixed_rim(vertices, faces, free2)
            vertices = smooth_region(vertices, faces, free2, edge_weights)

    return vertices, faces, patch_face_mask

smooth_region_boundary 

smooth_region_boundary(
    vertices: array[vec3],
    faces: array[int32],
    region: array[bool],
    iterations: int = 4,
    *,
    vertex_faces: tuple[array[int32], array[int32]]
    | None = None,
) -> wp.array[wp.vec3]

Straighten the rim curve of a face region, by sliding the vertices that lie on it.

The third member of the region family, and the one whose subject is the boundary itself. smooth_region smooths the surface across the rim and smooth_region_fixed_rim holds the rim still while smoothing inside it; both leave the rim's own path through the mesh exactly where the face selection put it. That path is usually ragged -- a selection by angle, by height or by a paint stroke follows triangle edges and zigzags -- and this is what makes it a smooth curve, without moving the surface off itself in any visible way.

The rim is smoothed as a level set rather than as a polyline. A field is pinned to -1 on the region's vertices and +1 outside, its harmonic interpolation is solved over the band of vertices that touch both sides, and each of those vertices is moved onto the zero level set of the result. Harmonic interpolation is smooth, so its zero set is a smooth curve; the vertices slide along the surface to sit on it. Repeating recomputes the field from the moved positions, and the move is damped so the sequence settles rather than oscillates.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
region array[bool]

(n_faces,) boolean mask selecting the region whose rim is smoothed. The rim is derived from it; only vertices that touch both a selected and an unselected face may move.

 required
iterations int

Number of solve-and-project passes. 0 returns a copy.

 4
vertex_faces tuple[array[int32], array[int32]] | None

Optional precomputed vertex_face_adjacency as (vertex_faces, offsets). Depends on the connectivity alone -- which this filter never changes -- so one CSR serves every call over the same mesh, and Trimesh.vertex_face_adjacency has it cached.

 None

Returns:

Type Description
array[vec3]

Positions on vertices.device with the rim band slid along the surface. Connectivity and the face mask are untouched, so region stays valid against the result.

Raises:

Type Description
ValueError

If iterations is negative, or if region is not a length-n_faces wp.bool array.
RuntimeError

If vertices, faces, region and vertex_faces are not all on one device.
See Also

smooth_region Smooths the surface across the rim, leaving the rim's path alone. smooth_region_fixed_rim Smooths inside the rim, holding it fixed. region_boundary_edges The rim as an edge list, which is what this leaves in a better place.

Notes

It slides vertices; it does not re-cut the rim

The rim can only become as smooth as the existing triangulation lets it: a vertex slides to the level set but the rim still passes through the same vertices, so a curve that wants to cross a triangle diagonally cannot. Retriangulating the band first (flip_by_objective, or split_faces_along_field to cut along the level set outright) is the way to get past that, and is deliberately not folded in here -- it changes the face buffer, which this promises not to.

Source code in triwarp/smoothing.py 
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
2030
2031
2032
2033
2034
2035
2036
2037
2038
2039
2040
2041
2042
2043
2044
2045
2046
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
2061
2062
2063
2064
2065
2066
2067
2068
2069
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079
2080
2081
2082
2083
2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108
2109
2110
2111
2112
2113
2114
2115
2116
2117
2118
def smooth_region_boundary(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    region: wp.array[wp.bool],
    iterations: int = 4,
    *,
    vertex_faces: tuple[wp.array[wp.int32], wp.array[wp.int32]] | None = None,
) -> wp.array[wp.vec3]:
    """
    Straighten the *rim curve* of a face region, by sliding the vertices that lie on it.

    The third member of the region family, and the one whose subject is the boundary itself.
    [`smooth_region`][triwarp.smoothing.smooth_region] smooths the surface *across* the rim and
    [`smooth_region_fixed_rim`][triwarp.smoothing.smooth_region_fixed_rim] holds the rim still while
    smoothing inside it; both leave the rim's own path through the mesh exactly where the face
    selection put it. That path is usually ragged -- a selection by angle, by height or by a paint
    stroke follows triangle edges and zigzags -- and this is what makes it a smooth curve, without
    moving the surface off itself in any visible way.

    The rim is smoothed as a **level set** rather than as a polyline. A field is pinned to ``-1`` on
    the region's vertices and ``+1`` outside, its harmonic interpolation is solved over the band of
    vertices that touch both sides, and each of those vertices is moved onto the zero level set of
    the result. Harmonic interpolation is smooth, so its zero set is a smooth curve; the vertices
    slide along the surface to sit on it. Repeating recomputes the field from the moved positions,
    and the move is damped so the sequence settles rather than oscillates.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    region
        ``(n_faces,)`` boolean mask selecting the region whose rim is smoothed. The rim is derived
        from it; only vertices that touch both a selected and an unselected face may move.
    iterations
        Number of solve-and-project passes. ``0`` returns a copy.
    vertex_faces
        Optional precomputed [`vertex_face_adjacency`][triwarp.adjacency.vertex_face_adjacency] as
        ``(vertex_faces, offsets)``. Depends on the connectivity alone -- which this filter never
        changes -- so one CSR serves every call over the same mesh, and
        [`Trimesh.vertex_face_adjacency`][triwarp.mesh.Trimesh.vertex_face_adjacency] has it cached.

    Returns
    -------
    wp.array[wp.vec3]
        Positions on ``vertices.device`` with the rim band slid along the surface. Connectivity and
        the face mask are untouched, so ``region`` stays valid against the result.

    Raises
    ------
    ValueError
        If ``iterations`` is negative, or if ``region`` is not a length-``n_faces`` ``wp.bool``
        array.
    RuntimeError
        If ``vertices``, ``faces``, ``region`` and ``vertex_faces`` are not all on one device.

    See Also
    --------
    [`smooth_region`][triwarp.smoothing.smooth_region]
        Smooths the surface across the rim, leaving the rim's path alone.
    [`smooth_region_fixed_rim`][triwarp.smoothing.smooth_region_fixed_rim]
        Smooths inside the rim, holding it fixed.
    [`region_boundary_edges`][triwarp.selection.region_boundary_edges]
        The rim as an edge list, which is what this leaves in a better place.

    Notes
    -----
    !!! note "It slides vertices; it does not re-cut the rim"
        The rim can only become as smooth as the *existing* triangulation lets it: a vertex slides
        to the level set but the rim still passes through the same vertices, so a curve that wants
        to cross a triangle diagonally cannot. Retriangulating the band first
        ([`flip_by_objective`][triwarp.remesh.flip_by_objective], or
        [`split_faces_along_field`][triwarp.intersection.split_faces_along_field] to cut along the
        level set outright) is the way to get past that, and is deliberately not folded in here --
        it changes the face buffer, which this promises not to.
    """
    require_same_device(vertices=vertices, faces=faces, region=region, vertex_faces=vertex_faces)
    if iterations < 0:
        raise ValueError(f"iterations must be non-negative, got {iterations}")
    device = faces.device
    n_faces = int(faces.shape[0]) // 3
    if len(region.shape) != 1 or region.shape[0] != n_faces or region.dtype is not wp.bool:
        raise ValueError(
            f"region must be a length-{n_faces} wp.bool array, got shape {tuple(region.shape)} "
            f"of {region.dtype}"
        )
    n_vertices = int(vertices.shape[0])
    if n_vertices == 0 or n_faces == 0 or iterations == 0:
        return wp.clone(vertices)

    inside = _incident_vertex_mask(faces, region, n_vertices)
    free = _region_rim_vertices(faces, region, inside)
    fixed_mask = wp.empty(n_vertices, dtype=wp.bool, device=device)
    wp.map(kernel_array.mask_not, free, out=fixed_mask)
    field = wp.empty((1, n_vertices), dtype=wp.float64, device=device)
    wp.map(kernel_smoothing.region_side_value, inside, out=field[0])

    vf_indices, offsets = (
        vertex_faces
        if vertex_faces is not None
        else tw.adjacency.vertex_face_adjacency(faces, n_vertices=n_vertices)
    )
    positions = wp.clone(vertices)
    nxt = wp.empty(n_vertices, dtype=wp.vec3, device=device)
    solved = wp.empty(n_vertices, dtype=wp.float64, device=device)
    for _ in range(iterations):
        # The cotangent weights depend on the positions, so the field is re-solved every pass; the
        # pinned +-1 entries are the same each time and only the band's values change.
        operator = wps.bsr_scale(laplacian.cotmatrix(positions, faces, dtype=wp.float64), -1.0)
        solution, free_map, n_free = twl.min_quad_with_fixed(
            operator, fixed_mask, twt.as_array2d(field, wp.float64)
        )
        if n_free == 0:
            break
        wp.copy(solved, field[0])
        wp.launch(
            kernel_smoothing.scatter_free_scalar,
            dim=n_vertices,
            inputs=[fixed_mask, free_map, solution[0], solved],
            device=device,
        )
        wp.launch(
            kernel_smoothing.project_to_zero_isoline,
            dim=n_vertices,
            inputs=[positions, faces, offsets, vf_indices, solved, free, _ISOLINE_DAMPING, nxt],
            device=device,
        )
        positions, nxt = nxt, positions
    return positions

filter_scalar_laplacian 

filter_scalar_laplacian(
    values: array[float32],
    vertices: array[vec3],
    faces: array[int32],
    lamb: float = 0.5,
    iterations: int = 10,
    laplacian_operator: BsrMatrix[float32] | None = None,
) -> wp.array[wp.float32]

Diffuse a per-vertex scalar field through the 1-ring averaging operator.

The scalar counterpart of filter_laplacian, running the identical explicit step q' = q + lamb (L q - q) on a field rather than on positions. Reach for it whenever a computed per-vertex quantity is noisier than the thing it will drive — curvature before it selects features, a heat or occlusion field before it becomes a weight, a per-vertex sizing function before it drives remeshing.

On a closed mesh, lamb=1.0, iterations=1 is exactly MeshLab's apply_scalar_smoothing_per_vertex: one unweighted 1-ring average with the vertex's own value excluded. The defaults instead match the rest of this module (a partial step, repeated), which is the gentler behaviour a caller usually wants.

It differs from MeshLab on a boundary

VCG's VertexQualityLaplacian smooths a boundary vertex along the boundary curve only — it averages just that vertex's two boundary neighbours and ignores the rest of its ring, so the boundary values evolve as an independent 1D field. That is a different operator, not this one with other constants, and it is why the pymeshlab oracle in tests/test_smoothing.py runs on closed fixtures. Here a boundary vertex averages its whole ring like any other; to hold the boundary fixed instead, restore its values after the call.

Parameters:

Name Type Description Default
values array[float32]

Length-n_vertices wp.float32 field to smooth. Not modified.

 required
vertices array[vec3]

(n_vertices,) mesh vertex positions — needed only to build the operator, and ignored when laplacian_operator is supplied.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
lamb float

Step size. 0 leaves the field unchanged, 1 replaces each value by its 1-ring average outright. Values above 1 overshoot and are unstable.

 0.5
iterations int

Number of passes. 0 returns a copy of the input.

 10
laplacian_operator BsrMatrix[float32] | None

Optional precomputed row-stochastic operator (see laplacian). Pass it to hoist the build out of a loop, or to smooth with inverse-edge-length weights instead.

When None it is built as laplacian(vertices, faces, symmetric=True) — the symmetric 1-ring, unlike filter_laplacian, which defaults to trimesh's directed mesh.edges adjacency. The two agree on a closed mesh and differ on a boundary vertex, where the directed adjacency is missing some of its neighbours; for a scalar field the symmetric ring is both the defensible choice and the one MeshLab makes.

 None

Returns:

Type Description
array[float32]

Length-n_vertices smoothed field on values.device.

Raises:

Type Description
ValueError

If values is not length n_vertices.
RuntimeError

If values, vertices and faces are not all on one device.
See Also

shortest_path_envelope filter_laplacian triwarp.laplacian.laplacian

Source code in triwarp/smoothing.py 
2180
2181
2182
2183
2184
2185
2186
2187
2188
2189
2190
2191
2192
2193
2194
2195
2196
2197
2198
2199
2200
2201
2202
2203
2204
2205
2206
2207
2208
2209
2210
2211
2212
2213
2214
2215
2216
2217
2218
2219
2220
2221
2222
2223
2224
2225
2226
2227
2228
2229
2230
2231
2232
2233
2234
2235
2236
2237
2238
2239
2240
2241
2242
2243
2244
2245
2246
2247
2248
2249
2250
2251
2252
2253
2254
2255
2256
2257
2258
2259
2260
2261
2262
2263
2264
2265
2266
2267
2268
2269
2270
2271
2272
2273
2274
2275
2276
2277
2278
2279
def filter_scalar_laplacian(
    values: wp.array[wp.float32],
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    lamb: float = 0.5,
    iterations: int = 10,
    laplacian_operator: wps.BsrMatrix[wp.float32] | None = None,
) -> wp.array[wp.float32]:
    """
    Diffuse a per-vertex scalar field through the 1-ring averaging operator.

    The scalar counterpart of [`filter_laplacian`][triwarp.smoothing.filter_laplacian], running the
    identical explicit step ``q' = q + lamb (L q - q)`` on a field rather than on positions. Reach
    for it whenever a computed per-vertex quantity is noisier than the thing it will drive —
    curvature before it selects features, a heat or occlusion field before it becomes a weight, a
    per-vertex sizing function before it drives remeshing.

    On a **closed** mesh, ``lamb=1.0, iterations=1`` is exactly MeshLab's
    ``apply_scalar_smoothing_per_vertex``: one unweighted 1-ring average with the vertex's own value
    excluded. The defaults instead match the rest of this module (a partial step, repeated), which
    is the gentler behaviour a caller usually wants.

    !!! note "It differs from MeshLab on a boundary"
        VCG's ``VertexQualityLaplacian`` smooths a boundary vertex along the **boundary curve
        only** — it averages just that vertex's two boundary neighbours and ignores the rest of its
        ring, so the boundary values evolve as an independent 1D field. That is a different
        operator, not this one with other constants, and it is why the pymeshlab oracle in
        ``tests/test_smoothing.py`` runs on closed fixtures. Here a boundary vertex averages its
        whole ring like any other; to hold the boundary fixed instead, restore its values after the
        call.

    Parameters
    ----------
    values
        Length-``n_vertices`` ``wp.float32`` field to smooth. Not modified.
    vertices
        ``(n_vertices,)`` mesh vertex positions — needed only to build the operator, and ignored
        when ``laplacian_operator`` is supplied.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    lamb
        Step size. ``0`` leaves the field unchanged, ``1`` replaces each value by its 1-ring
        average outright. Values above ``1`` overshoot and are unstable.
    iterations
        Number of passes. ``0`` returns a copy of the input.
    laplacian_operator
        Optional precomputed row-stochastic operator (see
        [`laplacian`][triwarp.laplacian.laplacian]). Pass it to hoist the build out of a loop, or
        to smooth with inverse-edge-length weights instead.

        When ``None`` it is built as ``laplacian(vertices, faces, symmetric=True)`` — the
        **symmetric** 1-ring, unlike [`filter_laplacian`][triwarp.smoothing.filter_laplacian],
        which defaults to trimesh's directed ``mesh.edges`` adjacency. The two agree on a closed
        mesh and differ on a boundary vertex, where the directed adjacency is missing some of its
        neighbours; for a scalar field the symmetric ring is both the defensible choice and the one
        MeshLab makes.

    Returns
    -------
    wp.array[wp.float32]
        Length-``n_vertices`` smoothed field on ``values.device``.

    Raises
    ------
    ValueError
        If ``values`` is not length ``n_vertices``.
    RuntimeError
        If ``values``, ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`shortest_path_envelope`][triwarp.graph.shortest_path_envelope]
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian]
    [`triwarp.laplacian.laplacian`][triwarp.laplacian.laplacian]
    """
    require_same_device(values=values, vertices=vertices, faces=faces)
    device = values.device
    n = int(vertices.shape[0])
    if int(values.shape[0]) != n:
        raise ValueError(
            f"values must have one entry per vertex, got {values.shape[0]} for {n} vertices"
        )

    out = wp.clone(values)
    if n == 0 or iterations == 0:
        return out

    operator = _resolved_operator(vertices, faces, laplacian_operator, symmetric=True)
    nxt = wp.empty(n, dtype=wp.float32, device=device)
    coeff = wp.float32(lamb)
    for _ in range(iterations):
        wp.launch(
            kernel_smoothing.diffuse_scalar_pass,
            dim=n,
            inputs=[operator.offsets, operator.columns, operator.values, out, coeff],
            outputs=[nxt],
            device=device,
        )
        out, nxt = nxt, out
    return out

filter_normals 

filter_normals(
    vertices: array[vec3],
    faces: array[int32],
    *,
    iterations: int = 20,
    threshold: float = 60.0,
    face_adjacency: Array2dInt32 | None = None,
    face_normals: array[vec3] | None = None,
    face_areas: array[float32] | None = None,
) -> wp.array[wp.vec3]

Smooth the face normals without moving a vertex, preserving creases.

Each pass replaces a face's normal with the area-weighted average of its own and those of its edge-neighbours whose normal is within threshold of it. That gate is what makes this feature-preserving rather than isotropic: across a crease the two normals disagree by more than the threshold and never average, so a sharp edge survives any number of passes while noise on a flat region diffuses away in a few. MeshLab's apply_normal_smoothing_per_face.

The result is a normal field that is no longer the geometric normal of any triangle — it is the target for filter_two_step's second half, and useful on its own for shading.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions. Not modified.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
iterations int

Number of diffusion passes over the normal field. MeshLab's stepnormalnum, whose default of 20 is this one.

 20
threshold float

Angle in degrees beyond which two neighbouring faces refuse to average. 0 averages nothing and 180 averages everything (isotropic). MeshLab's normalthr, default 60. Must be in [0, 180].

 60.0
face_adjacency Array2dInt32 | None

Optional (m, 2) adjacency from face_adjacency; recomputed when None. Pass it to hoist the build out of a loop — the topology never changes here, so it is safe to reuse.

 None
face_normals array[vec3] | None

Optional length-n_faces unit face normals and matching areas from face_normals_and_areas; recomputed together when either is None. Trimesh.face_normals and Trimesh.face_areas cache the pair.

 None
face_areas array[float32] | None

See face_normals.

 None

Returns:

Type Description
array[vec3]

Length-n_faces unit normals on vertices.device. A fresh buffer even when face_normals is passed: the passes rewrite the field in place.

Raises:

Type Description
ValueError

If threshold is outside [0, 180].
RuntimeError

If vertices, faces, face_adjacency, face_normals and face_areas are not all on one device.
See Also

filter_two_step triwarp.triangles.face_normals_and_areas

Source code in triwarp/smoothing.py 
2282
2283
2284
2285
2286
2287
2288
2289
2290
2291
2292
2293
2294
2295
2296
2297
2298
2299
2300
2301
2302
2303
2304
2305
2306
2307
2308
2309
2310
2311
2312
2313
2314
2315
2316
2317
2318
2319
2320
2321
2322
2323
2324
2325
2326
2327
2328
2329
2330
2331
2332
2333
2334
2335
2336
2337
2338
2339
2340
2341
2342
2343
2344
2345
2346
2347
2348
2349
2350
2351
2352
2353
2354
2355
2356
2357
2358
2359
2360
2361
2362
2363
2364
2365
2366
2367
2368
2369
2370
2371
2372
2373
2374
2375
2376
2377
2378
2379
2380
2381
2382
2383
2384
2385
2386
2387
2388
2389
2390
2391
2392
2393
2394
2395
2396
2397
2398
2399
2400
2401
2402
2403
2404
2405
2406
def filter_normals(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    *,
    iterations: int = 20,
    threshold: float = 60.0,
    face_adjacency: twt.Array2dInt32 | None = None,
    face_normals: wp.array[wp.vec3] | None = None,
    face_areas: wp.array[wp.float32] | None = None,
) -> wp.array[wp.vec3]:
    """
    Smooth the *face normals* without moving a vertex, preserving creases.

    Each pass replaces a face's normal with the area-weighted average of its own and those of its
    edge-neighbours **whose normal is within ``threshold`` of it**. That gate is what makes this
    feature-preserving rather than isotropic: across a crease the two normals disagree by more than
    the threshold and never average, so a sharp edge survives any number of passes while noise on a
    flat region diffuses away in a few. MeshLab's ``apply_normal_smoothing_per_face``.

    The result is a normal field that is no longer the geometric normal of any triangle — it is the
    *target* for [`filter_two_step`][triwarp.smoothing.filter_two_step]'s second half, and useful on
    its own for shading.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions. Not modified.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    iterations
        Number of diffusion passes over the normal field. MeshLab's ``stepnormalnum``, whose default
        of ``20`` is this one.
    threshold
        Angle in **degrees** beyond which two neighbouring faces refuse to average. ``0`` averages
        nothing and ``180`` averages everything (isotropic). MeshLab's ``normalthr``, default
        ``60``. Must be in ``[0, 180]``.
    face_adjacency
        Optional ``(m, 2)`` adjacency from
        [`face_adjacency`][triwarp.adjacency.face_adjacency]; recomputed when ``None``. Pass it to
        hoist the build out of a loop — the topology never changes here, so it is safe to reuse.
    face_normals
        Optional length-``n_faces`` unit face normals and matching areas from
        [`face_normals_and_areas`][triwarp.triangles.face_normals_and_areas]; recomputed together
        when either is ``None``. [`Trimesh.face_normals`][triwarp.mesh.Trimesh.face_normals] and
        [`Trimesh.face_areas`][triwarp.mesh.Trimesh.face_areas] cache the pair.
    face_areas
        See ``face_normals``.

    Returns
    -------
    wp.array[wp.vec3]
        Length-``n_faces`` unit normals on ``vertices.device``. A fresh buffer even when
        ``face_normals`` is passed: the passes rewrite the field in place.

    Raises
    ------
    ValueError
        If ``threshold`` is outside ``[0, 180]``.
    RuntimeError
        If ``vertices``, ``faces``, ``face_adjacency``, ``face_normals`` and ``face_areas`` are not
        all on one device.

    See Also
    --------
    [`filter_two_step`][triwarp.smoothing.filter_two_step]
    [`triwarp.triangles.face_normals_and_areas`][triwarp.triangles.face_normals_and_areas]
    """
    require_same_device(
        vertices=vertices,
        faces=faces,
        face_adjacency=face_adjacency,
        face_normals=face_normals,
        face_areas=face_areas,
    )
    if not 0.0 <= threshold <= 180.0:
        raise ValueError(f"threshold must be in [0, 180] degrees, got {threshold}")

    device = vertices.device
    n_faces = int(faces.shape[0]) // 3
    if face_normals is None or face_areas is None:
        face_normals, face_areas = face_normals_and_areas(vertices, faces)
    # Each pass rewrites the normal field in place, so a caller-supplied buffer is copied rather
    # than clobbered -- ``Trimesh.face_normals`` is one array shared with every other reader.
    normals, areas = wp.clone(face_normals), face_areas
    if n_faces == 0 or iterations <= 0:
        return normals

    if face_adjacency is None:
        face_adjacency = tw.adjacency.face_adjacency(faces, n_vertices=int(vertices.shape[0]))
    m = int(face_adjacency.shape[0])
    threshold_cos = wp.float32(math.cos(math.radians(threshold)))

    accumulated = wp.empty(n_faces, dtype=wp.vec3, device=device)
    # The map is hoisted out of the pass loop: a cached ``wp.map`` call re-resolves its kernel in
    # Python every time, which adds up over many passes.
    seed = wp.map(
        kernel_smoothing.seed_weighted_normal, normals, areas, out=accumulated, return_kernel=True
    )
    renormalize = wp.map(wp.normalize, accumulated, out=normals, return_kernel=True)
    # Two launches per pass rather than three. A pass is seed -> crease-gated neighbour
    # accumulation -> normalize, and the scatter in the middle needs the whole seeded buffer, so the
    # fusable pair is the normalization with the *next* pass's seed. Peeling the first seed off the
    # front is what puts them next to each other.
    wp.launch(seed, dim=n_faces, inputs=[normals, areas], outputs=[accumulated], device=device)
    for index in range(iterations):
        if m > 0:
            wp.launch(
                kernel_smoothing.accumulate_smoothed_normals,
                dim=m,
                inputs=[normals, areas, face_adjacency, threshold_cos, accumulated],
                device=device,
            )
        if index == iterations - 1:
            wp.launch(
                renormalize, dim=n_faces, inputs=[accumulated], outputs=[normals], device=device
            )
        else:
            wp.launch(
                kernel_smoothing.renormalize_and_reseed,
                dim=n_faces,
                inputs=[areas, accumulated],
                outputs=[normals],
                device=device,
            )
    return normals

filter_two_step 

filter_two_step(
    vertices: array[vec3],
    faces: array[int32],
    *,
    iterations: int = 3,
    threshold: float = 60.0,
    normal_iterations: int = 20,
    fit_iterations: int = 20,
) -> wp.array[wp.vec3]

Feature-preserving smoothing: filter the face normals, then fit the vertices to them.

The one filter in this module that does not blur a crease. Every other scheme here diffuses positions directly, which cannot distinguish noise from a sharp edge — both are high-frequency. Two-step smoothing (Ohtake et al., and MeshLab's apply_coord_two_steps_smoothing) separates the two questions:

  1. Where should the surface face? Diffuse the normal field with filter_normals, whose threshold refuses to average across a crease. Noise is removed; the crease is not.
  2. Where should the vertices go? Move each vertex so the planes through its incident faces' centroids with those filtered normals agree as well as possible — a gradient step repeated fit_iterations times, with no free step size (the average over incident faces is the step).

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions. Not modified.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer. The topology never changes.

 required
iterations int

Number of outer passes; each one re-derives the normals from the current positions and refits. MeshLab's stepsmoothnum, default 3.

 3
threshold float

Crease threshold in degrees for the normal filter. See filter_normals. MeshLab's normalthr, default 60.

 60.0
normal_iterations int

Normal-diffusion passes per outer pass. MeshLab's stepnormalnum, default 20.

 20
fit_iterations int

Vertex-fitting gradient steps per outer pass. MeshLab's stepfitnum, default 20.

 20

Returns:

Type Description
array[vec3]

Smoothed (n_vertices,) vertex positions on vertices.device.

Raises:

Type Description
ValueError

If threshold is outside [0, 180].
RuntimeError

If vertices and faces are not all on one device.
See Also

filter_normals filter_laplacian filter_sharpen

Notes

There is no volume constraint and none is needed: the fitting step moves each vertex along the filtered normals rather than toward its neighbours' mean, so the systematic inward drift that filter_laplacian has to correct for does not arise. A flat region is already a fixed point of both halves.

Source code in triwarp/smoothing.py 
2409
2410
2411
2412
2413
2414
2415
2416
2417
2418
2419
2420
2421
2422
2423
2424
2425
2426
2427
2428
2429
2430
2431
2432
2433
2434
2435
2436
2437
2438
2439
2440
2441
2442
2443
2444
2445
2446
2447
2448
2449
2450
2451
2452
2453
2454
2455
2456
2457
2458
2459
2460
2461
2462
2463
2464
2465
2466
2467
2468
2469
2470
2471
2472
2473
2474
2475
2476
2477
2478
2479
2480
2481
2482
2483
2484
2485
2486
2487
2488
2489
2490
2491
2492
2493
2494
2495
2496
2497
2498
2499
2500
2501
2502
2503
2504
2505
2506
2507
2508
def filter_two_step(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    *,
    iterations: int = 3,
    threshold: float = 60.0,
    normal_iterations: int = 20,
    fit_iterations: int = 20,
) -> wp.array[wp.vec3]:
    """
    Feature-preserving smoothing: filter the face normals, then fit the vertices to them.

    The one filter in this module that does **not** blur a crease. Every other scheme here diffuses
    positions directly, which cannot distinguish noise from a sharp edge — both are high-frequency.
    Two-step smoothing (Ohtake et al., and MeshLab's ``apply_coord_two_steps_smoothing``) separates
    the two questions:

    1. **Where should the surface face?** Diffuse the *normal* field with
       [`filter_normals`][triwarp.smoothing.filter_normals], whose threshold refuses to average
       across a crease. Noise is removed; the crease is not.
    2. **Where should the vertices go?** Move each vertex so the planes through its incident faces'
       centroids with those filtered normals agree as well as possible — a gradient step repeated
       ``fit_iterations`` times, with no free step size (the average over incident faces is the
       step).

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions. Not modified.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer. The topology never changes.
    iterations
        Number of outer passes; each one re-derives the normals from the current positions and
        refits. MeshLab's ``stepsmoothnum``, default ``3``.
    threshold
        Crease threshold in **degrees** for the normal filter. See
        [`filter_normals`][triwarp.smoothing.filter_normals]. MeshLab's ``normalthr``, default
        ``60``.
    normal_iterations
        Normal-diffusion passes per outer pass. MeshLab's ``stepnormalnum``, default ``20``.
    fit_iterations
        Vertex-fitting gradient steps per outer pass. MeshLab's ``stepfitnum``, default ``20``.

    Returns
    -------
    wp.array[wp.vec3]
        Smoothed ``(n_vertices,)`` vertex positions on ``vertices.device``.

    Raises
    ------
    ValueError
        If ``threshold`` is outside ``[0, 180]``.
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`filter_normals`][triwarp.smoothing.filter_normals]
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian]
    [`filter_sharpen`][triwarp.smoothing.filter_sharpen]

    Notes
    -----
    There is no volume constraint and none is needed: the fitting step moves each vertex *along* the
    filtered normals rather than toward its neighbours' mean, so the systematic inward drift that
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian] has to correct for does not arise. A
    flat region is already a fixed point of both halves.
    """
    require_same_device(vertices=vertices, faces=faces)
    device = vertices.device
    n = int(vertices.shape[0])
    n_faces = int(faces.shape[0]) // 3
    out = wp.clone(vertices)
    if n == 0 or n_faces == 0 or iterations <= 0:
        return out

    # The topology is fixed, so the adjacency is built once for every pass of both halves.
    adjacency = tw.adjacency.face_adjacency(faces, n_vertices=n)
    delta = wp.empty(n, dtype=wp.vec3, device=device)
    counts = wp.empty(n, dtype=wp.float32, device=device)
    # Hoisted out of the doubly-nested pass loop, where it would otherwise run
    # ``iterations * fit_iterations`` times.
    fit_step = wp.map(
        kernel_smoothing.apply_fit_step, out, delta, counts, out=out, return_kernel=True
    )
    for _ in range(iterations):
        normals = filter_normals(
            out, faces, iterations=normal_iterations, threshold=threshold, face_adjacency=adjacency
        )
        for _fit in range(fit_iterations):
            delta.zero_()
            counts.zero_()
            wp.launch(
                kernel_smoothing.fit_vertices_to_normals,
                dim=n_faces,
                inputs=[out, faces, normals, delta, counts],
                device=device,
            )
            wp.launch(fit_step, dim=n, inputs=[out, delta, counts], outputs=[out], device=device)
    return out

filter_sharpen 

filter_sharpen(
    vertices: array[vec3],
    faces: array[int32],
    *,
    weight: float = 0.3,
    weight_original: float = 1.0,
    iterations: int = 5,
    laplacian_operator: BsrMatrix[float32] | None = None,
) -> wp.array[wp.vec3]

Sharpen a surface by adding back the detail a smoothing pass removes.

The inverse of a smoothing filter, built from one: the difference between the mesh and its Laplacian-smoothed self is its high-frequency content, so adding a multiple of that difference exaggerates every feature. This is unsharp masking, the photographic technique, applied to vertex positions instead of pixels — MeshLab spells it apply_coord_unsharp_mask — and it is the standard way to recover crispness lost to an earlier smoothing or a decimation.

The name says what it returns, which is positions rather than a mask: the *_mask family in this package is boolean throughout.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions. Not modified.

 required
faces array[int32]

Length-3 * n_faces wp.int32 triangle index buffer.

 required
weight float

How much of the detail to add back. 0 leaves the mesh unchanged and large values amplify noise as readily as features. MeshLab's weight, default 0.3.

 0.3
weight_original float

Multiplier on the original position, MeshLab's weightorig. 1 (the default) keeps the surface where it is and only sharpens; anything else scales the whole mesh about the origin.

 1.0
iterations int

Laplacian passes used to define "smoothed", so this sets the scale of detail that gets amplified: more passes remove lower frequencies and therefore sharpen more coarsely. MeshLab's iterations, default 5.

 5
laplacian_operator BsrMatrix[float32] | None

Optional precomputed operator, forwarded to filter_laplacian.

 None

Returns:

Type Description
array[vec3]

Sharpened (n_vertices,) vertex positions on vertices.device.

Raises:

Type Description
RuntimeError

If vertices and faces are not all on one device.
See Also

filter_laplacian filter_two_step

Source code in triwarp/smoothing.py 
2511
2512
2513
2514
2515
2516
2517
2518
2519
2520
2521
2522
2523
2524
2525
2526
2527
2528
2529
2530
2531
2532
2533
2534
2535
2536
2537
2538
2539
2540
2541
2542
2543
2544
2545
2546
2547
2548
2549
2550
2551
2552
2553
2554
2555
2556
2557
2558
2559
2560
2561
2562
2563
2564
2565
2566
2567
2568
2569
2570
2571
2572
2573
2574
2575
2576
2577
2578
2579
2580
2581
2582
2583
2584
2585
2586
2587
2588
2589
2590
2591
def filter_sharpen(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    *,
    weight: float = 0.3,
    weight_original: float = 1.0,
    iterations: int = 5,
    laplacian_operator: wps.BsrMatrix[wp.float32] | None = None,
) -> wp.array[wp.vec3]:
    """
    Sharpen a surface by adding back the detail a smoothing pass removes.

    The inverse of a smoothing filter, built from one: the difference between the mesh and its
    Laplacian-smoothed self *is* its high-frequency content, so adding a multiple of that difference
    exaggerates every feature. This is **unsharp masking**, the photographic technique, applied to
    vertex positions instead of pixels — MeshLab spells it ``apply_coord_unsharp_mask`` — and it is
    the standard way to recover crispness lost to an earlier smoothing or a decimation.

    The name says what it returns, which is *positions* rather than a mask: the ``*_mask`` family
    in this package is boolean throughout.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions. Not modified.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    weight
        How much of the detail to add back. ``0`` leaves the mesh unchanged and large values amplify
        noise as readily as features. MeshLab's ``weight``, default ``0.3``.
    weight_original
        Multiplier on the original position, MeshLab's ``weightorig``. ``1`` (the default) keeps the
        surface where it is and only sharpens; anything else scales the whole mesh about the origin.
    iterations
        Laplacian passes used to define "smoothed", so this sets the *scale* of detail that gets
        amplified: more passes remove lower frequencies and therefore sharpen more coarsely.
        MeshLab's ``iterations``, default ``5``.
    laplacian_operator
        Optional precomputed operator, forwarded to
        [`filter_laplacian`][triwarp.smoothing.filter_laplacian].

    Returns
    -------
    wp.array[wp.vec3]
        Sharpened ``(n_vertices,)`` vertex positions on ``vertices.device``.

    Raises
    ------
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian]
    [`filter_two_step`][triwarp.smoothing.filter_two_step]
    """
    require_same_device(vertices=vertices, faces=faces)
    device = vertices.device
    n = int(vertices.shape[0])
    out = wp.empty(n, dtype=wp.vec3, device=device)
    if n == 0:
        return out
    # No volume constraint on the inner smoothing: its rescaling would leak into the difference and
    # show up as a uniform scaling of the sharpened mesh rather than as detail.
    smoothed = filter_laplacian(
        vertices,
        faces,
        lamb=1.0,
        iterations=iterations,
        volume_constraint=False,
        laplacian_operator=laplacian_operator,
    )
    wp.map(
        kernel_smoothing.unsharp_step,
        vertices,
        smoothed,
        wp.float32(weight),
        wp.float32(weight_original),
        out=out,
    )
    return out