Skip to content

geodesic_walk

triwarp.geodesic_walk

Straightest geodesics: walk a fixed distance across a mesh in a fixed direction.

A straightest geodesic is what you get by walking forward and, at every edge, unfolding the two incident triangles into a common plane and continuing in a straight line. It is the surface analogue of "go that way for this far", which makes it the tool for exponential maps, streamline tracing and extending a direction field along a surface — and unlike a shortest path it is fixed by an initial condition rather than by two endpoints.

Both entry points are batched over many rays: one thread walks one ray, and the traced polylines come back packed into one buffer with CSR offsets, the same shape boundary_loops_batched uses.

trace_from_vertex is the exponential map of the surface, and its inverse is log_map: one takes a tangent direction and a distance to a point, the other takes a point back to the direction and distance that reach it. The two live apart because they are different methods -- this one unfolds triangles combinatorially, that one solves a vector-heat system -- and geodesic distance by the heat method is a third, heat_geodesic.

Cone points

A path that runs exactly into a vertex has no unique straightest continuation — the angle around a vertex is not 2 * pi, so "straight through" is ambiguous. Such a crossing is resolved by unfolding across the edge the walk reached the vertex along: the correct limit for a path passing arbitrarily close to the vertex, but not the split-the-angle convention. Paths through high-curvature vertices therefore drift from geometry-central's by about the angle defect.

trace_from_vertex 

trace_from_vertex(
    vertices: array[vec3],
    faces: array[int32],
    start_vertices: array[int32],
    directions: array[vec3],
    twins: array[int32] | None = None,
    rings: tuple[array[int32], array[int32], array[bool]]
    | None = None,
    frames: tuple[array[vec3], array[vec3], array[vec3]]
    | None = None,
    max_steps: int = _DEFAULT_MAX_STEPS,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

Trace a straightest geodesic from each of a batch of vertices.

Each ray starts at vertices[start_vertices[r]] and walks in directions[r], for an arc length equal to that direction's length after projection into the vertex's tangent plane — a unit direction traces at most a unit distance, and a direction along the normal traces nothing. That is potpourri3d.GeodesicTracer's convention too.

Which incident face the ray starts in is decided in the vertex's flattened tangent space (halfedge_tangent_angles): rescaling the incident corner angles to a full turn makes the fan a disk, so every tangent direction lands in exactly one wedge, including directions a naive per-face projection would place outside all of them. At a boundary vertex the fan spans only half a disk, and a direction outside it — pointing off the surface — traces nothing.

A walk stops early when it reaches the mesh boundary or exceeds max_steps edge crossings.

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
start_vertices array[int32]

(n_rays,) wp.int32 start vertex per ray.

 required
directions array[vec3]

(n_rays,) initial directions; the length of each sets how far its ray is traced.

 required
twins array[int32] | None

Optional precomputed halfedge_twins.

 None
rings tuple[array[int32], array[int32], array[bool]] | None

Optional precomputed vertex_one_rings, used to find each ray's starting face.

 None
frames tuple[array[vec3], array[vec3], array[vec3]] | None

Optional precomputed vertex_tangent_frames. They define each vertex's tangent plane, which sets both the trace length and the starting wedge.

 None
max_steps int

Maximum edge crossings per ray.

 _DEFAULT_MAX_STEPS

Returns:

Name Type Description
points array[vec3]

All traced points, packed ray after ray, on vertices.device.
offsets array[int32]

Length n_rays + 1; ray r owns points[offsets[r] : offsets[r + 1]], beginning at its start vertex. A ray always contributes at least one point.

Raises:

Type Description
RuntimeError

If vertices, faces, start_vertices, directions, twins, rings and frames are not all on one device.
See Also

trace_from_face trace_polylines heat_geodesic

Source code in triwarp/geodesic_walk.py 
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 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
def trace_from_vertex(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    start_vertices: wp.array[wp.int32],
    directions: wp.array[wp.vec3],
    twins: wp.array[wp.int32] | None = None,
    rings: tuple[wp.array[wp.int32], wp.array[wp.int32], wp.array[wp.bool]] | None = None,
    frames: tuple[wp.array[wp.vec3], wp.array[wp.vec3], wp.array[wp.vec3]] | None = None,
    max_steps: int = _DEFAULT_MAX_STEPS,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Trace a straightest geodesic from each of a batch of vertices.

    Each ray starts at ``vertices[start_vertices[r]]`` and walks in ``directions[r]``, for an arc
    length equal to that direction's length *after projection into the vertex's tangent plane* — a
    unit direction traces at most a unit distance, and a direction along the normal traces nothing.
    That is ``potpourri3d.GeodesicTracer``'s convention too.

    Which incident face the ray starts in is decided in the vertex's *flattened* tangent space
    ([`halfedge_tangent_angles`][triwarp.tangent_space.halfedge_tangent_angles]): rescaling the
    incident corner angles to a full turn makes the fan a disk, so every tangent direction lands in
    exactly one wedge, including directions a naive per-face projection would place outside all of
    them. At a boundary vertex the fan spans only half a disk, and a direction outside it — pointing
    off the surface — traces nothing.

    A walk stops early when it reaches the mesh boundary or exceeds ``max_steps`` edge crossings.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    start_vertices
        ``(n_rays,)`` ``wp.int32`` start vertex per ray.
    directions
        ``(n_rays,)`` initial directions; the length of each sets how far its ray is traced.
    twins
        Optional precomputed [`halfedge_twins`][triwarp.halfedge.halfedge_twins].
    rings
        Optional precomputed [`vertex_one_rings`][triwarp.halfedge.vertex_one_rings], used to find
        each ray's starting face.
    frames
        Optional precomputed [`vertex_tangent_frames`][triwarp.tangent_space.vertex_tangent_frames].
        They define each vertex's tangent plane, which sets both the trace length and the starting
        wedge.
    max_steps
        Maximum edge crossings per ray.

    Returns
    -------
    points : wp.array[wp.vec3]
        All traced points, packed ray after ray, on ``vertices.device``.
    offsets : wp.array[wp.int32]
        Length ``n_rays + 1``; ray ``r`` owns ``points[offsets[r] : offsets[r + 1]]``, beginning at
        its start vertex. A ray always contributes at least one point.

    Raises
    ------
    RuntimeError
        If ``vertices``, ``faces``, ``start_vertices``, ``directions``, ``twins``, ``rings`` and
        ``frames`` are not all on one device.

    See Also
    --------
    [`trace_from_face`][triwarp.geodesic_walk.trace_from_face]
    [`trace_polylines`][triwarp.geodesic_walk.trace_polylines]
    [`heat_geodesic`][triwarp.heat.heat_geodesic]
    """
    require_same_device(
        vertices=vertices,
        faces=faces,
        start_vertices=start_vertices,
        directions=directions,
        twins=twins,
        rings=rings,
        frames=frames,
    )
    device = vertices.device
    n_rays = int(start_vertices.shape[0])
    n_vertices = int(vertices.shape[0])
    if n_rays == 0:
        return wp.empty(0, dtype=wp.vec3, device=device), wp.zeros(
            max(n_rays + 1, 1), dtype=wp.int32, device=device
        )
    # No short-circuit on ``faces.shape[0] == 0``: unlike ``trace_from_face``, a start *vertex* is
    # still a valid reference into ``vertices`` with no faces at all (every vertex is isolated), and
    # the general path below already produces the documented one-point-per-ray answer for that case
    # -- ``vertex_one_rings`` gives every vertex an empty ring, and the kernel's
    # ``start_direction_at_vertex`` then reports "isolated vertex, nowhere to go" (``f == -1``),
    # which ``trace_from_vertices`` turns into exactly one written point, the start vertex itself.
    # A short-circuit here returned an empty polyline per ray instead, contradicting this function's
    # own "a ray always contributes at least one point" guarantee.

    if twins is None:
        twins = halfedge_twins(faces, n_vertices=n_vertices)
    if rings is None:
        rings = vertex_one_rings(faces, twins=twins, n_vertices=n_vertices)
    ring_halfedges, ring_offsets, is_boundary = rings
    if frames is None:
        frames = tw.tangent_space.vertex_tangent_frames(vertices, faces, rings=rings)
    basis_x, basis_y, normals = frames

    inputs = [
        vertices,
        faces,
        twins,
        tw.triangles.face_angles(vertices, faces),
        ring_offsets,
        ring_halfedges,
        is_boundary,
        basis_x,
        basis_y,
        normals,
        start_vertices,
        directions,
        wp.int32(max_steps),
        wp.float32(_length_epsilon(vertices, faces)),
    ]
    return _trace(kernel_geodesic_walk.trace_from_vertices, inputs, n_rays, device)

trace_from_face 

trace_from_face(
    vertices: array[vec3],
    faces: array[int32],
    start_faces: array[int32],
    start_barycentric: array[vec3],
    directions: array[vec3],
    twins: array[int32] | None = None,
    max_steps: int = _DEFAULT_MAX_STEPS,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

Trace a straightest geodesic from each of a batch of barycentric points inside faces.

As trace_from_vertex, but each ray starts at an interior point of a known face, so no wedge search is needed. This is the form for streamlines of a face-based vector field.

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
start_faces array[int32]

(n_rays,) wp.int32 starting face per ray.

 required
start_barycentric array[vec3]

(n_rays,) barycentric coordinates of each start point within its face.

 required
directions array[vec3]

(n_rays,) initial directions; the length of each sets how far its ray is traced.

 required
twins array[int32] | None

Optional precomputed halfedge_twins.

 None
max_steps int

Maximum edge crossings per ray.

 _DEFAULT_MAX_STEPS

Returns:

Name Type Description
points array[vec3]

All traced points, packed ray after ray, on vertices.device.
offsets array[int32]

Length n_rays + 1 CSR bounds into points.

Raises:

Type Description
RuntimeError

If vertices, faces, start_faces, start_barycentric, directions and twins are not all on one device.
See Also

trace_from_vertex barycentric_to_points

Source code in triwarp/geodesic_walk.py 
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
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
def trace_from_face(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    start_faces: wp.array[wp.int32],
    start_barycentric: wp.array[wp.vec3],
    directions: wp.array[wp.vec3],
    twins: wp.array[wp.int32] | None = None,
    max_steps: int = _DEFAULT_MAX_STEPS,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Trace a straightest geodesic from each of a batch of barycentric points inside faces.

    As [`trace_from_vertex`][triwarp.geodesic_walk.trace_from_vertex], but each ray
    starts at an interior point of a known face, so no wedge search is needed. This is the form for
    streamlines of a face-based vector field.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    start_faces
        ``(n_rays,)`` ``wp.int32`` starting face per ray.
    start_barycentric
        ``(n_rays,)`` barycentric coordinates of each start point within its face.
    directions
        ``(n_rays,)`` initial directions; the length of each sets how far its ray is traced.
    twins
        Optional precomputed [`halfedge_twins`][triwarp.halfedge.halfedge_twins].
    max_steps
        Maximum edge crossings per ray.

    Returns
    -------
    points : wp.array[wp.vec3]
        All traced points, packed ray after ray, on ``vertices.device``.
    offsets : wp.array[wp.int32]
        Length ``n_rays + 1`` CSR bounds into ``points``.

    Raises
    ------
    RuntimeError
        If ``vertices``, ``faces``, ``start_faces``, ``start_barycentric``, ``directions`` and
        ``twins`` are not all on one device.

    See Also
    --------
    [`trace_from_vertex`][triwarp.geodesic_walk.trace_from_vertex]
    [`barycentric_to_points`][triwarp.triangles.barycentric_to_points]
    """
    require_same_device(
        vertices=vertices,
        faces=faces,
        start_faces=start_faces,
        start_barycentric=start_barycentric,
        directions=directions,
        twins=twins,
    )
    device = vertices.device
    n_rays = int(start_faces.shape[0])
    if n_rays == 0 or int(faces.shape[0]) == 0:
        return wp.empty(0, dtype=wp.vec3, device=device), wp.zeros(
            max(n_rays + 1, 1), dtype=wp.int32, device=device
        )

    if twins is None:
        twins = halfedge_twins(faces, n_vertices=int(vertices.shape[0]))

    inputs = [
        vertices,
        faces,
        twins,
        start_faces,
        start_barycentric,
        directions,
        wp.int32(max_steps),
        wp.float32(_length_epsilon(vertices, faces)),
    ]
    return _trace(kernel_geodesic_walk.trace_from_faces, inputs, n_rays, device)

descend_field 

descend_field(
    vertices: array[vec3],
    faces: array[int32],
    values: array[float64],
    starts: array[int32],
    *,
    stop_value: float = 0.0,
    twins: array[int32] | None = None,
    vertex_faces: tuple[array[int32], array[int32]]
    | None = None,
    gradients: array[vec3d] | None = None,
    max_steps: int = _DEFAULT_MAX_STEPS,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

Follow a per-vertex scalar field downhill from each of a batch of start vertices.

The field-descent counterpart of this module's two straightest-geodesic tracers: where those are fixed by an initial direction, this one is fixed by a field, and it goes wherever that field decreases fastest. Fed a geodesic distance field it traces the geodesic back to its source, which is what geodesic_path is; fed any other scalar it traces that scalar's flow lines.

Three cases, and the field value at each written point strictly decreases in all of them -- which is what makes the walk terminate rather than orbit:

  1. Inside a face the piecewise-linear interpolant's gradient is constant, so the path is a straight segment to the exit edge.
  2. Along an edge, when the face across it has a descent that points back: the walk slides to the edge's lower-valued endpoint.
  3. At a vertex, where the field has no single gradient: each incident face is asked whether its own descent direction points into the fan wedge, and the steepest admissible one wins.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions.

 required
faces array[int32]

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

 required
values array[float64]

Length-n_vertices wp.float64 field to descend. float64 because the fields this serves are exponentially decaying -- see triwarp.laplacian.face_gradients, which computes the per-face gradient this walks along.

 required
starts array[int32]

(n_paths,) wp.int32 vertices to descend from, one path each.

 required
stop_value float

Field value at which a path stops. The default 0.0 is what a distance field's source sits at.

 0.0
twins array[int32] | None

Optional precomputed triwarp.halfedge.halfedge_twins.

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

Optional precomputed triwarp.adjacency.vertex_face_adjacency pair, which case 3 needs. (vertex_faces, offsets), values first, as that function returns it and as every packed pair in the package is spelled -- passing it the other way round reads offsets as face indices and raises nothing, since both are wp.int32.

 None
gradients array[vec3d] | None

Optional precomputed per-face gradient of values. Pass it when descending the same field from several batches.

 None
max_steps int

Cap on steps per path. A path that hits it is returned truncated rather than reported.

 _DEFAULT_MAX_STEPS

Returns:

Type Description
(points, offsets)

points holds every path's polyline end to end and offsets is the length-n_paths + 1 CSR bound, the same packing trace_from_vertex returns and trace_polylines slices.

Raises:

Type Description
ValueError

If values does not have one entry per vertex, or gradients is given and does not have one entry per face.
RuntimeError

If vertices, faces, values, starts, twins, vertex_faces and gradients are not all on one device.
Notes

A path can stop before reaching stop_value, and the caller can tell: its last point is not within tolerance of a vertex whose value is at the stop. That happens at a local minimum of the field, on a flat face where the gradient vanishes, at the mesh boundary, and when max_steps runs out. None of those is an error -- a field with several minima has several basins, and this walks the one it starts in.

See Also

geodesic_path The distance-field case, which is what this is usually reached for. trace_from_vertex The direction-driven walk, for a straightest geodesic rather than a shortest one. triwarp.laplacian.face_gradients

Source code in triwarp/geodesic_walk.py 
250
251
252
253
254
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
def descend_field(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    values: wp.array[wp.float64],
    starts: wp.array[wp.int32],
    *,
    stop_value: float = 0.0,
    twins: wp.array[wp.int32] | None = None,
    vertex_faces: tuple[wp.array[wp.int32], wp.array[wp.int32]] | None = None,
    gradients: wp.array[wp.vec3d] | None = None,
    max_steps: int = _DEFAULT_MAX_STEPS,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Follow a per-vertex scalar field downhill from each of a batch of start vertices.

    The **field-descent** counterpart of this module's two straightest-geodesic tracers: where those
    are fixed by an initial direction, this one is fixed by a *field*, and it goes wherever that
    field decreases fastest. Fed a geodesic distance field it traces the geodesic back to its
    source, which is what [`geodesic_path`][triwarp.geodesic_walk.geodesic_path] is; fed any other
    scalar it traces that scalar's flow lines.

    Three cases, and the field value at each written point strictly decreases in all of them --
    which is what makes the walk terminate rather than orbit:

    1. **Inside a face** the piecewise-linear interpolant's gradient is constant, so the path is a
       straight segment to the exit edge.
    2. **Along an edge**, when the face across it has a descent that points back: the walk slides to
       the edge's lower-valued endpoint.
    3. **At a vertex**, where the field has no single gradient: each incident face is asked whether
       its own descent direction points into the fan wedge, and the steepest admissible one wins.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer.
    values
        Length-``n_vertices`` ``wp.float64`` field to descend. ``float64`` because the fields this
        serves are exponentially decaying -- see
        [`triwarp.laplacian.face_gradients`][triwarp.laplacian.face_gradients], which computes the
        per-face gradient this walks along.
    starts
        ``(n_paths,)`` ``wp.int32`` vertices to descend from, one path each.
    stop_value
        Field value at which a path stops. The default ``0.0`` is what a distance field's source
        sits at.
    twins
        Optional precomputed [`triwarp.halfedge.halfedge_twins`][triwarp.halfedge.halfedge_twins].
    vertex_faces
        Optional precomputed
        [`triwarp.adjacency.vertex_face_adjacency`][triwarp.adjacency.vertex_face_adjacency] pair,
        which case 3 needs. ``(vertex_faces, offsets)``, values first, as that function returns it
        and as every packed pair in the package is spelled -- passing it the other way round reads
        offsets as face indices and raises nothing, since both are ``wp.int32``.
    gradients
        Optional precomputed per-face gradient of ``values``. Pass it when descending the same field
        from several batches.
    max_steps
        Cap on steps per path. A path that hits it is returned truncated rather than reported.

    Returns
    -------
    points, offsets
        ``points`` holds every path's polyline end to end and ``offsets`` is the
        length-``n_paths + 1`` CSR bound, the same packing
        [`trace_from_vertex`][triwarp.geodesic_walk.trace_from_vertex] returns and
        [`trace_polylines`][triwarp.geodesic_walk.trace_polylines] slices.

    Raises
    ------
    ValueError
        If ``values`` does not have one entry per vertex, or ``gradients`` is given and does not
        have one entry per face.
    RuntimeError
        If ``vertices``, ``faces``, ``values``, ``starts``, ``twins``, ``vertex_faces`` and
        ``gradients`` are not all on one device.

    Notes
    -----
    A path can stop before reaching ``stop_value``, and the caller can tell: its last point is not
    within tolerance of a vertex whose value is at the stop. That happens at a **local minimum** of
    the field, on a flat face where the gradient vanishes, at the mesh **boundary**, and when
    ``max_steps`` runs out. None of those is an error -- a field with several minima has several
    basins, and this walks the one it starts in.

    See Also
    --------
    [`geodesic_path`][triwarp.geodesic_walk.geodesic_path]
        The distance-field case, which is what this is usually reached for.
    [`trace_from_vertex`][triwarp.geodesic_walk.trace_from_vertex]
        The direction-driven walk, for a *straightest* geodesic rather than a shortest one.
    [`triwarp.laplacian.face_gradients`][triwarp.laplacian.face_gradients]
    """
    require_same_device(
        vertices=vertices,
        faces=faces,
        values=values,
        starts=starts,
        twins=twins,
        vertex_faces=vertex_faces,
        gradients=gradients,
    )
    device = vertices.device
    n_vertices = int(vertices.shape[0])
    if int(values.shape[0]) != n_vertices:
        raise ValueError(
            f"values must have one entry per vertex, got {values.shape[0]} for {n_vertices}"
        )
    n_faces = int(faces.shape[0]) // 3
    if gradients is not None and int(gradients.shape[0]) != n_faces:
        # Unlike ``values``, a mismatched ``gradients`` was an out-of-bounds read with no
        # exception -- ``descent_walk``/``descend_at_vertex`` index it as ``gradients[f]`` for
        # every face index up to ``n_faces - 1`` (silent host-heap corruption on CPU, garbage or a
        # crash on CUDA), the exact hazard a caller-supplied precomputed cache exists to avoid.
        raise ValueError(
            f"gradients must have one entry per face, got {gradients.shape[0]} for {n_faces}"
        )
    n_paths = int(starts.shape[0])
    if n_paths == 0:
        return (
            wp.empty(0, dtype=wp.vec3, device=device),
            wp.zeros(1, dtype=wp.int32, device=device),
        )

    if twins is None:
        twins = halfedge_twins(faces, n_vertices=n_vertices)
    if vertex_faces is None:
        vertex_faces = tw.adjacency.vertex_face_adjacency(faces, n_vertices=n_vertices)
    if gradients is None:
        gradients = tw.laplacian.face_gradients(vertices, faces, values)
    incident_faces, face_offsets = vertex_faces

    inputs = [
        vertices,
        faces,
        twins,
        face_offsets,
        incident_faces,
        values,
        gradients,
        starts,
        wp.float64(stop_value),
        wp.int32(max_steps),
        wp.float32(_length_epsilon(vertices, faces)),
    ]
    return _trace(kernel_geodesic_walk.descent_paths, inputs, n_paths, device)

geodesic_path 

geodesic_path(
    vertices: array[vec3],
    faces: array[int32],
    source: array[int32],
    targets: array[int32],
    *,
    t: float | None = None,
    operators: object | None = None,
    max_steps: int = _DEFAULT_MAX_STEPS,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]

Trace a path across the surface from each target vertex back to a source set.

The point-to-point geodesic, as against this module's straightest walks and heat_geodesic's distance field: one heat solve gives the distance to the source everywhere, and descending it from a target follows that geodesic back. Every target shares the one solve, so a thousand paths to one source cost one system and a thousand independent walks -- which is why the signature is one source and many targets rather than a list of pairs.

Parameters:

Name Type Description Default
vertices array[vec3]

(n_vertices,) mesh vertex positions.

 required
faces array[int32]

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

 required
source array[int32]

(k,) wp.int32 source vertices. Several make the paths run to whichever is nearest, since the field is the distance to the set.

 required
targets array[int32]

(n_paths,) wp.int32 vertices to trace from.

 required
t float | None

Heat diffusion time, forwarded to heat_geodesic. None uses its default.

 None
operators object | None

Prebuilt HeatOperators for this mesh, to spare the factorization when several sources are traced on one mesh.

 None
max_steps int

Cap on steps per path.

 _DEFAULT_MAX_STEPS

Returns:

Type Description
(points, offsets)

Packed polylines and their CSR bounds, each running from its target to the source. Slice with trace_polylines and measure with triwarp.polyline.polyline_length.

Raises:

Type Description
RuntimeError

If vertices, faces, source and targets are not all on one device.

Examples:

source = tw.array.arange(1, device=v.device)
targets = tw.array.arange(int(v.shape[0]), device=v.device)
points, offsets = tw.geodesic_walk.geodesic_path(v, f, source, targets)
Notes

The path is as accurate as the field it descends, and no more. The heat method's distance is first-order, so this is an approximate geodesic: measured against potpourri3d's edge-flip geodesics -- which are exact -- the length comes out a few per cent long, and the excess is the field's error rather than the walk's. It is never shorter than the true geodesic, which is the invariant worth testing against.

A path that cannot reach the source stops early rather than failing; see descend_field for the four ways that happens.

See Also

descend_field The walk itself, for descending any other scalar field. triwarp.heat.heat_geodesic The field, when the distance is wanted and not the path.

Source code in triwarp/geodesic_walk.py 
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
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
def geodesic_path(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    source: wp.array[wp.int32],
    targets: wp.array[wp.int32],
    *,
    t: float | None = None,
    operators: object | None = None,
    max_steps: int = _DEFAULT_MAX_STEPS,
) -> tuple[wp.array[wp.vec3], wp.array[wp.int32]]:
    """
    Trace a path across the surface from each target vertex back to a source set.

    The **point-to-point** geodesic, as against this module's straightest walks and
    [`heat_geodesic`][triwarp.heat.heat_geodesic]'s distance *field*: one heat solve gives
    the distance to the source everywhere, and descending it from a target follows that geodesic
    back. Every target shares the one solve, so a thousand paths to one source cost one system and a
    thousand independent walks -- which is why the signature is one source and many targets rather
    than a list of pairs.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` flat triangle index buffer.
    source
        ``(k,)`` ``wp.int32`` source vertices. Several make the paths run to whichever is nearest,
        since the field is the distance to the *set*.
    targets
        ``(n_paths,)`` ``wp.int32`` vertices to trace from.
    t
        Heat diffusion time, forwarded to
        [`heat_geodesic`][triwarp.heat.heat_geodesic]. ``None`` uses its default.
    operators
        Prebuilt [`HeatOperators`][triwarp.heat.HeatOperators] for this mesh, to spare the
        factorization when several sources are traced on one mesh.
    max_steps
        Cap on steps per path.

    Returns
    -------
    points, offsets
        Packed polylines and their CSR bounds, each running **from its target to the source**. Slice
        with [`trace_polylines`][triwarp.geodesic_walk.trace_polylines] and measure with
        [`triwarp.polyline.polyline_length`][triwarp.polyline.polyline_length].

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

    Examples
    --------
    ```python
    source = tw.array.arange(1, device=v.device)
    targets = tw.array.arange(int(v.shape[0]), device=v.device)
    points, offsets = tw.geodesic_walk.geodesic_path(v, f, source, targets)
    ```

    Notes
    -----
    **The path is as accurate as the field it descends, and no more.** The heat method's distance is
    first-order, so this is an *approximate* geodesic: measured against ``potpourri3d``'s edge-flip
    geodesics -- which are exact -- the length comes out a few per cent long, and the excess is the
    field's error rather than the walk's. It is never *shorter* than the true geodesic, which is the
    invariant worth testing against.

    A path that cannot reach the source stops early rather than failing; see
    [`descend_field`][triwarp.geodesic_walk.descend_field] for the four ways that happens.

    See Also
    --------
    [`descend_field`][triwarp.geodesic_walk.descend_field]
        The walk itself, for descending any other scalar field.
    [`triwarp.heat.heat_geodesic`][triwarp.heat.heat_geodesic]
        The field, when the distance is wanted and not the path.
    """
    require_same_device(
        vertices=vertices, faces=faces, source=source, targets=targets, operators=operators
    )
    distance = tw.heat.heat_geodesic(vertices, faces, source, t, operators)  # type: ignore[arg-type]
    return descend_field(vertices, faces, distance, targets, stop_value=0.0, max_steps=max_steps)

shorten_loop 

shorten_loop(
    vertices: array[vec3],
    faces: array[int32],
    loops: Sequence[array[int32]],
    *,
    max_iter: int = 100,
    tolerance: float = 0.0,
    twins: array[int32] | None = None,
    rings: tuple[array[int32], array[int32], array[bool]]
    | None = None,
) -> tuple[list[wp.array[wp.int32]], int]

Shorten closed edge loops within their homotopy class, keeping them on mesh edges.

Takes exactly what homology_generators returns -- a list of vertex-index cycles -- and returns cycles of the same kind, shorter. The loops a tree-cotree construction produces are as long and as jagged as the spanning trees that built them, which is fine for a basis and useless as a curve; this makes them short enough to look at, cut along, or measure.

Each sweep rewrites the loop locally. Around one of its vertices b, with neighbours a and c on the loop, the sub-path a -> b -> c is replaced by whichever way round b's link is shorter, when either beats going through b. Both alternatives lie in the star of b, which is a disk, so the replacement cannot change the loop's homotopy class -- and since it is only ever accepted when it is strictly shorter, the total length falls monotonically. A loop that doubles back on itself has the spur contracted away by the same rule.

Shorter, not geodesic

The result stays on the edge graph, so it is a local minimum over edge paths rather than the shortest curve on the surface -- reaching that means letting the loop cross face interiors, which needs an intrinsic triangulation it can flip. The gap to the true geodesic is small, and is widest where the mesh is a regular grid whose rows are not geodesics, because no one-ring move can step the loop off a row without lengthening the edge path first -- so the result is a true local minimum over edge paths, and the remaining gap is the edge graph's, not the sweep's.

Sweeps ratchet the loop one position at a time, so the count needed to converge grows with the loop's own length, which is what max_iter has to cover.

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
loops Sequence[array[int32]]

Closed vertex-index cycles, each a wp.int32 array whose consecutive entries share an edge, as do its last and first. Loops shorter than three vertices are returned unchanged.

 required
max_iter int

Cap on the number of sweeps. Two sweeps of opposite parity are needed to give every position a turn, so an odd cap leaves one parity class one visit short. Sweeping stops as soon as two consecutive sweeps accept nothing, which is what makes the default generous rather than expensive.

 100
tolerance float

Absolute length a replacement must save to be accepted. The default 0.0 accepts any strict improvement, which is what makes the result independent of the sweep count; raise it to stop the last few sweeps chasing float32 noise on a fine mesh.

 0.0
twins array[int32] | None

Optional precomputed halfedge_twins.

 None
rings tuple[array[int32], array[int32], array[bool]] | None

Optional precomputed vertex_one_rings as (ring_halfedges, offsets, is_boundary). Depends on the connectivity alone, so one CSR serves every fan walk over the same mesh -- Trimesh.vertex_one_rings has it cached, and passing it skips the vertex-manifold check and the host readback that check costs.

 None

Returns:

Name Type Description
loops list[array[int32]]

One shortened cycle per input loop, in the same order, on faces.device.
sweeps int

How many sweeps ran. Below max_iter this means the loops stopped changing, so the answer is locally minimal; equal to it, the cap bound the result.

Raises:

Type Description
TypeError

If any loop is not a rank-1 wp.int32 array.
RuntimeError

If vertices, faces, loops, twins and rings are not all on one device.
See Also

homology_generators Produces the loops this shortens. polyline_length Measures the result, after gathering the positions. geodesic_path The open, endpoint-to-endpoint problem, solved by descending a heat field instead.

Source code in triwarp/geodesic_walk.py 
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
565
566
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
def shorten_loop(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    loops: Sequence[wp.array[wp.int32]],
    *,
    max_iter: int = 100,
    tolerance: float = 0.0,
    twins: wp.array[wp.int32] | None = None,
    rings: tuple[wp.array[wp.int32], wp.array[wp.int32], wp.array[wp.bool]] | None = None,
) -> tuple[list[wp.array[wp.int32]], int]:
    """
    Shorten closed edge loops within their homotopy class, keeping them on mesh edges.

    Takes exactly what [`homology_generators`][triwarp.homology.homology_generators] returns -- a
    list of vertex-index cycles -- and returns cycles of the same kind, shorter. The loops a
    tree-cotree construction produces are as long and as jagged as the spanning trees that built
    them, which is fine for a *basis* and useless as a curve; this makes them short enough to look
    at, cut along, or measure.

    Each sweep rewrites the loop *locally*. Around one of its vertices ``b``, with neighbours ``a``
    and ``c`` on the loop, the sub-path ``a -> b -> c`` is replaced by whichever way round ``b``'s
    link is shorter, when either beats going through ``b``. Both alternatives lie in the star of
    ``b``, which is a disk, so the replacement cannot change the loop's homotopy class -- and since
    it is only ever accepted when it is strictly shorter, the total length falls monotonically. A
    loop that doubles back on itself has the spur contracted away by the same rule.

    !!! note "Shorter, not geodesic"
        The result stays **on the edge graph**, so it is a local minimum over edge paths rather than
        the shortest curve on the surface -- reaching that means letting the loop cross face
        interiors, which needs an intrinsic triangulation it can flip. The gap to the true geodesic
        is small, and is widest where the mesh is a regular grid whose rows are not geodesics,
        because no one-ring move can step the loop off a row without lengthening the edge path
        first -- so the result is a true local minimum over edge paths, and the remaining gap is the
        edge graph's, not the sweep's.

        Sweeps ratchet the loop one position at a time, so the count needed to converge grows with
        the loop's own length, which is what ``max_iter`` has to cover.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    loops
        Closed vertex-index cycles, each a ``wp.int32`` array whose consecutive entries share an
        edge, as do its last and first. Loops shorter than three vertices are returned unchanged.
    max_iter
        Cap on the number of sweeps. Two sweeps of opposite parity are needed to give every
        position a turn, so an odd cap leaves one parity class one visit short. Sweeping stops as
        soon as two consecutive sweeps accept nothing, which is what makes the default generous
        rather than expensive.
    tolerance
        Absolute length a replacement must save to be accepted. The default ``0.0`` accepts any
        strict improvement, which is what makes the result independent of the sweep count; raise it
        to stop the last few sweeps chasing float32 noise on a fine mesh.
    twins
        Optional precomputed [`halfedge_twins`][triwarp.halfedge.halfedge_twins].
    rings
        Optional precomputed [`vertex_one_rings`][triwarp.halfedge.vertex_one_rings] as
        ``(ring_halfedges, offsets, is_boundary)``. Depends on the connectivity alone, so one CSR
        serves every fan walk over the same mesh --
        [`Trimesh.vertex_one_rings`][triwarp.mesh.Trimesh.vertex_one_rings] has it cached, and
        passing it skips the vertex-manifold check and the host readback that check costs.

    Returns
    -------
    loops : list[wp.array[wp.int32]]
        One shortened cycle per input loop, in the same order, on ``faces.device``.
    sweeps : int
        How many sweeps ran. Below ``max_iter`` this means the loops stopped changing, so the
        answer is locally minimal; equal to it, the cap bound the result.

    Raises
    ------
    TypeError
        If any loop is not a rank-1 ``wp.int32`` array.
    RuntimeError
        If ``vertices``, ``faces``, ``loops``, ``twins`` and ``rings`` are not all on one device.

    See Also
    --------
    [`homology_generators`][triwarp.homology.homology_generators]
        Produces the loops this shortens.
    [`polyline_length`][triwarp.polyline.polyline_length]
        Measures the result, after gathering the positions.
    [`geodesic_path`][triwarp.geodesic_walk.geodesic_path]
        The open, endpoint-to-endpoint problem, solved by descending a heat field instead.
    """
    require_same_device(vertices=vertices, faces=faces, loops=loops, twins=twins, rings=rings)
    device = faces.device
    loops = list(loops)
    for loop in loops:
        twt.ensure_ndim(loop, 1, dtype=wp.int32)
    if not loops or int(faces.shape[0]) == 0 or max_iter <= 0:
        return loops, 0

    n_vertices = int(vertices.shape[0])
    if twins is None:
        twins = halfedge_twins(faces, n_vertices=n_vertices)
    ring_halfedges, ring_offsets, is_boundary = (
        rings if rings is not None else vertex_one_rings(faces, twins=twins, n_vertices=n_vertices)
    )

    # ``copy=False``: each sweep reads ``packed`` and writes a freshly sized buffer, so the
    # first pack can alias the caller's loops. A sweep that accepts nothing leaves them alone.
    packed, starts = tw.array.pack_1d_arrays(loops, copy=False)
    # The offsets `pack_1d_arrays` returns are not total-terminated, and every kernel below reads
    # `loop_offsets[l + 1]`, so terminate them once here rather than special-casing the last loop.
    loop_offsets = tw.array.concatenate(
        [starts, wp.array([packed.shape[0]], dtype=wp.int32, device=device)]
    )
    n_loops = len(loops)
    changed = wp.zeros(1, dtype=wp.int32, device=device)

    sweeps = 0
    # Requires TWO CONSECUTIVE sweeps to accept nothing, not "this sweep is odd and accepted
    # nothing" -- a rewrite can shift which loop positions land on even/odd parity, so a single
    # unchanged sweep says nothing about whether the *other* parity would still find something.
    consecutive_unchanged = 0
    for sweep in range(max_iter):
        n_positions = int(packed.shape[0])
        if n_positions == 0:
            break
        position_loop = wp.empty(n_positions, dtype=wp.int32, device=device)
        wp.launch(
            kernel_array.segment_owner_labels,
            dim=n_loops,
            inputs=[loop_offsets, position_loop],
            device=device,
        )
        counts = wp.empty(n_positions, dtype=wp.int32, device=device)
        arc_slot = wp.empty(n_positions, dtype=wp.int32, device=device)
        arc_step = wp.empty(n_positions, dtype=wp.int32, device=device)
        changed.zero_()
        wp.launch(
            kernel_geodesic_walk.shorten_loop_counts,
            dim=n_positions,
            inputs=[
                vertices,
                faces,
                ring_offsets,
                ring_halfedges,
                is_boundary,
                packed,
                position_loop,
                loop_offsets,
                sweep % 2,
                tolerance,
                counts,
                arc_slot,
                arc_step,
                changed,
            ],
            device=device,
        )
        sweeps = sweep + 1
        # One 4-byte readback per sweep, and the only way to stop early: whether any replacement was
        # accepted is a device-side fact, and the alternative -- always running `max_iter` sweeps --
        # costs a full pass over every loop for each one that would have been skipped.
        if int(read_scalar(changed, 0)) == 0:
            consecutive_unchanged += 1
            if consecutive_unchanged >= 2:
                break  # both parities have now had a turn with nothing to do
            continue
        consecutive_unchanged = 0
        packed, loop_offsets = _rewrite_loops(
            faces, ring_offsets, ring_halfedges, packed, counts, arc_slot, arc_step, loop_offsets
        )
        packed, loop_offsets = _compact_repeats(packed, loop_offsets, n_loops)

    return tw.array.split(packed, loop_offsets[:n_loops]), sweeps

trace_polylines 

trace_polylines(
    points: array[vec3],
    offsets: array[int32],
    *,
    copy: bool = False,
) -> list[wp.array[wp.vec3]]

Slice a packed trace result into one polyline per ray.

The returned arrays are views

Each polyline slices the packed buffer, so holding one keeps them all alive and writing into one writes into the shared allocation. Pass copy=True for independent buffers.

Parameters:

Name Type Description Default
points array[vec3]

Packed traced points from trace_from_vertex or trace_from_face.

 required
offsets array[int32]

The matching length-n_rays + 1 CSR bounds.

 required
copy bool

Return independent buffers instead of views.

 False

Returns:

Type Description
list[array[vec3]]

One open polyline per ray, in ray order.

Raises:

Type Description
RuntimeError

If points and offsets are not all on one device.
See Also

polyline_length boundary_loops

Source code in triwarp/geodesic_walk.py 
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
787
788
789
790
791
def trace_polylines(
    points: wp.array[wp.vec3], offsets: wp.array[wp.int32], *, copy: bool = False
) -> list[wp.array[wp.vec3]]:
    """
    Slice a packed trace result into one polyline per ray.

    !!! note "The returned arrays are views"
        Each polyline slices the packed buffer, so holding one keeps them all alive and writing into
        one writes into the shared allocation. Pass ``copy=True`` for independent buffers.

    Parameters
    ----------
    points
        Packed traced points from
        [`trace_from_vertex`][triwarp.geodesic_walk.trace_from_vertex] or
        [`trace_from_face`][triwarp.geodesic_walk.trace_from_face].
    offsets
        The matching length-``n_rays + 1`` CSR bounds.
    copy
        Return independent buffers instead of views.

    Returns
    -------
    list[wp.array[wp.vec3]]
        One open polyline per ray, in ray order.

    Raises
    ------
    RuntimeError
        If ``points`` and ``offsets`` are not all on one device.

    See Also
    --------
    [`polyline_length`][triwarp.polyline.polyline_length]
    [`boundary_loops`][triwarp.boundary.boundary_loops]
    """
    require_same_device(points=points, offsets=offsets)
    # ``offsets`` is the total-terminated n + 1 form, and ``split`` wants the length-n one, whose
    # last segment already runs to the end of ``points``. The guard is required rather than
    # defensive: a no-ray trace returns a length-1 offsets array, and Warp rejects the resulting
    # zero-length slice outright ("Invalid indexing in slice: 0:0:1").
    if int(offsets.shape[0]) <= 1:
        return []
    return tw.array.split(points, offsets[:-1], copy=copy)