From a8d033e4a29c1d6a3b56efe6d8ec481ee999638c Mon Sep 17 00:00:00 2001 From: Andras Lasso Date: Wed, 9 Sep 2026 23:19:22 -0400 Subject: [PATCH 1/2] ENH: Support variable stent radius profiles in SDFStent deployment Generalize the capsule-chain stent SDF from a single shared radius (pill shape) to an arbitrary radius profile along the centerline: capsule_sdf, smin_sdf_capsule_contact_sculpt and compute_sdf_contact_displacements now accept the stent radius either as a scalar (unchanged behavior) or as a per-vertex array of shape (V,), linearly interpolated along each capsule segment (a chain of tapered capsules). This enables flared (funnel or trumpet shaped) ends, tapers, and other variable-radius stents. The per-segment distance/direction computation shared by the two SDF functions is factored into a common helper. A new cap_height_fraction parameter (default 1.0 = the original spherical caps) flattens the capsule end caps into half ellipsoids whose axial semi-axis is the given fraction of the local radius, implemented by scaling up the axial overshoot in each segment's closest-point computation. Flattened caps let a radius profile express concave features - a wide capsule's spherical end cap would otherwise bulge a full radius deep into a neighboring narrow region and wash the narrowing out - and keep a flared end from inflating a large ball into the vessel beyond the stent end. New geometry helpers build per-vertex radius profiles for a stent axis: stent_radius_profile (from arbitrary (position, radius) control points) and flared_stent_radius_profile (smoothstep flare at one end), both exported from svmorph.core. stent_bounding_box accepts a per-vertex radius array (its maximum is used for culling). All defaults preserve the previous behavior exactly; existing call sites are unchanged. Adds unit tests for radius interpolation, scalar/array equivalence, cap shapes, concave (dumbbell) profile preservation, and the profile helpers. Co-Authored-By: Claude Fable 5 --- svmorph/core/__init__.py | 6 +- svmorph/core/deformation.py | 145 +++++++++++++++++++++++++++--------- svmorph/core/geometry.py | 87 ++++++++++++++++++++++ tests/test_deformation.py | 92 +++++++++++++++++++++++ tests/test_geometry.py | 44 ++++++++++- 5 files changed, 334 insertions(+), 40 deletions(-) create mode 100644 tests/test_deformation.py diff --git a/svmorph/core/__init__.py b/svmorph/core/__init__.py index 7b042dbc..0a07cf84 100644 --- a/svmorph/core/__init__.py +++ b/svmorph/core/__init__.py @@ -11,7 +11,11 @@ set_node_indices, stent_bounding_box, ) -from svmorph.core.geometry import resample_stent_axis +from svmorph.core.geometry import ( + flared_stent_radius_profile, + resample_stent_axis, + stent_radius_profile, +) from svmorph.core.mesh_data import ( apply_displacements, compute_material_constants, diff --git a/svmorph/core/deformation.py b/svmorph/core/deformation.py index 271d0252..098f2a27 100644 --- a/svmorph/core/deformation.py +++ b/svmorph/core/deformation.py @@ -174,12 +174,13 @@ def compute_min_dist_and_direction(d: jx.Array, direction: jx.Array) -> tuple[jx return final_d, final_dir @jx.jit -def capsule_sdf(p: jx.Array, stent_vertices: jx.Array, r: float) -> jx.Array: - """Evaluate the signed distance field of a capsule-chain stent. +def capsule_sdf(p: jx.Array, stent_vertices: jx.Array, r: float | jx.Array, + cap_height_fraction: float = 1.0) -> jx.Array: + """Evaluate the signed distance field of a (tapered) capsule-chain stent. - Each consecutive pair of *stent_vertices* defines a capsule segment - with radius *r*. The SDF is reduced via smooth-minimum so the - iso-surface is C¹-continuous at segment junctions. + Each consecutive pair of *stent_vertices* defines a capsule segment. + The SDF is reduced via smooth-minimum so the iso-surface is + C¹-continuous at segment junctions. Parameters ---------- @@ -187,27 +188,83 @@ def capsule_sdf(p: jx.Array, stent_vertices: jx.Array, r: float) -> jx.Array: Query points, shape ``(N, 3)``. stent_vertices : jx.Array Stent axis vertices, shape ``(V, 3)``. - r : float - Capsule radius. + r : float or jx.Array + Capsule radius: scalar, or per-vertex radii of shape ``(V,)``, + linearly interpolated along each segment (arbitrary radius profile + along the centerline, e.g. a flared end). + cap_height_fraction : float + Axial semi-axis of the capsule end caps as a fraction of the local + radius; ``1.0`` (default) gives spherical caps, smaller values + flatten them into half ellipsoids (see + :func:`_tapered_capsule_segment_distances`). Returns ------- jx.Array Signed distance for each query point, shape ``(N, 1)``. """ + dist_to_surface_all, direction_all = _tapered_capsule_segment_distances( + p, stent_vertices, r, cap_height_fraction) + final_dist_to_surface, _ = jx.vmap(compute_min_dist_and_direction)(dist_to_surface_all, direction_all) + final_dist_to_surface = final_dist_to_surface[:, None] + return final_dist_to_surface + +def _tapered_capsule_segment_distances( + p: jx.Array, stent_vertices: jx.Array, r: float | jx.Array, cap_height_fraction: float, +) -> tuple[jx.Array, jx.Array]: + """Per-segment signed distances and outward directions of a tapered capsule chain. + + Each consecutive pair of *stent_vertices* defines a capsule segment. The + radius may be a single scalar (uniform stent) or a per-vertex array of + shape ``(V,)`` that is linearly interpolated along each segment, so the + stent can follow an arbitrary radius profile along the centerline (e.g. a + flared/trumpet end, a taper, or a local narrowing). + + *cap_height_fraction* controls the shape of the capsule end caps: the + axial overshoot beyond the segment ends (nonzero only where the + closest-point parameter clamps) is scaled up by its reciprocal, which + turns each spherical cap into a half ellipsoid whose axial semi-axis is + ``cap_height_fraction`` times the local radius. ``1.0`` gives the plain + spherical caps. Flattened caps let a radius profile express concave + features — a wide capsule's spherical cap would otherwise bulge a full + radius deep into a neighboring narrow region — and keep a flared stent + end from inflating a ball into the vessel beyond the stent end. + + Parameters + ---------- + p : jx.Array + Query points, shape ``(N, 3)``. + stent_vertices : jx.Array + Stent axis vertices, shape ``(V, 3)``. + r : float or jx.Array + Capsule radius: scalar, or per-vertex radii of shape ``(V,)``. + cap_height_fraction : float + Axial semi-axis of the capsule end caps as a fraction of the local + radius, in ``(0, 1]``. + + Returns + ------- + dist_to_surface_all : jx.Array + Per-segment signed distances, shape ``(N, num_segments)``. + direction_all : jx.Array + Per-segment outward direction vectors, shape ``(N, num_segments, 3)``. + """ + r = jnp.broadcast_to(jnp.asarray(r), (stent_vertices.shape[0],)) ba_all = jnp.diff(stent_vertices, axis=0) pa_all = p[:, None, :] - stent_vertices[None, :-1, :] ba_dot_pa_all = jnp.sum(pa_all * ba_all[None, :, :], axis=-1) ba_dot_ba_all = jnp.sum(ba_all**2, axis=-1) - h_all = jnp.clip(ba_dot_pa_all / ba_dot_ba_all, 0, 1) - axis_to_point_all = pa_all - h_all[:, :, None] * ba_all[None, :, :] + h_unclamped_all = ba_dot_pa_all / ba_dot_ba_all + h_all = jnp.clip(h_unclamped_all, 0, 1) + cap_axial_scale = 1.0 / cap_height_fraction + axis_to_point_all = (pa_all - h_all[:, :, None] * ba_all[None, :, :] + + ((h_unclamped_all - h_all) * (cap_axial_scale - 1.0))[:, :, None] * ba_all[None, :, :]) dist_all = jnp.linalg.norm(axis_to_point_all, axis=-1)[..., None] direction_all = axis_to_point_all / dist_all dist_all_squeezed = jnp.squeeze(dist_all, axis=-1) # shape: (num_mesh_points, num_segments) - dist_to_surface_all = dist_all_squeezed - r - final_dist_to_surface, _ = jx.vmap(compute_min_dist_and_direction)(dist_to_surface_all, direction_all) - final_dist_to_surface = final_dist_to_surface[:, None] - return final_dist_to_surface + r_at_closest_point_all = r[None, :-1] + h_all * (r[1:] - r[:-1])[None, :] + dist_to_surface_all = dist_all_squeezed - r_at_closest_point_all + return dist_to_surface_all, direction_all def kelvinlets_truncated_spherical_contraction( rv: jx.Array, f_scale: float, @@ -285,7 +342,8 @@ def kelvinlets_truncated_spherical_expansion( @jx.jit def smin_sdf_capsule_contact_sculpt( rv: jx.Array, stent_vertices: jx.Array, - r_current: float, + r_current: float | jx.Array, + cap_height_fraction: float = 1.0, ) -> tuple[jx.Array, jx.Array]: """Compute smooth-min SDF distances and outward directions from a capsule-chain stent. @@ -300,8 +358,15 @@ def smin_sdf_capsule_contact_sculpt( Query points, shape ``(N, 3)``. stent_vertices : jx.Array Stent axis vertices, shape ``(V, 3)``. - r_current : float - Current stent deployment radius. + r_current : float or jx.Array + Current stent deployment radius: scalar, or per-vertex radii of + shape ``(V,)``, linearly interpolated along each segment (arbitrary + radius profile along the centerline, e.g. a flared end). + cap_height_fraction : float + Axial semi-axis of the capsule end caps as a fraction of the local + radius; ``1.0`` (default) gives spherical caps, smaller values + flatten them into half ellipsoids (see + :func:`_tapered_capsule_segment_distances`). Returns ------- @@ -310,16 +375,8 @@ def smin_sdf_capsule_contact_sculpt( final_direction : jx.Array Unit outward direction from the stent axis, shape ``(N, 3)``. """ - ba_all = jnp.diff(stent_vertices, axis=0) - pa_all = rv[:, None, :] - stent_vertices[None, :-1, :] - ba_dot_pa_all = jnp.sum(pa_all * ba_all[None, :, :], axis=-1) - ba_dot_ba_all = jnp.sum(ba_all**2, axis=-1) - h_all = jnp.clip(ba_dot_pa_all / ba_dot_ba_all, 0, 1) - axis_to_point_all = pa_all - h_all[:, :, None] * ba_all[None, :, :] - dist_all = jnp.linalg.norm(axis_to_point_all, axis=-1)[..., None] - direction_all = axis_to_point_all / dist_all - dist_all_squeezed = jnp.squeeze(dist_all, axis=-1) # shape: (num_mesh_points, num_segments) - dist_to_surface_all = dist_all_squeezed - r_current + dist_to_surface_all, direction_all = _tapered_capsule_segment_distances( + rv, stent_vertices, r_current, cap_height_fraction) # Vectorize the folding over all mesh points: final_dist_to_surface, final_direction = jx.vmap(compute_min_dist_and_direction)(dist_to_surface_all, direction_all) final_dist_to_surface = final_dist_to_surface[:, None] @@ -560,8 +617,9 @@ def stent_bounding_box( Mesh vertices, shape ``(N, 3)``. stent_vertices : jx.Array Stent axis vertices, shape ``(V, 3)``. - target_stent_radius : float - Target stent radius. + target_stent_radius : float or jx.Array + Target stent radius (the maximum is used when a per-vertex radius + array is given). influence_radius : float Additional radial padding for the influence zone. contact_distance : float @@ -572,17 +630,18 @@ def stent_bounding_box( jx.Array Boolean mask of length *N*. """ - min_coords = jnp.min(stent_vertices, axis=0) - target_stent_radius - influence_radius - contact_distance - 0.01 * L() - max_coords = jnp.max(stent_vertices, axis=0) + target_stent_radius + influence_radius + contact_distance + 0.01 * L() + max_stent_radius = jnp.max(jnp.asarray(target_stent_radius)) + min_coords = jnp.min(stent_vertices, axis=0) - max_stent_radius - influence_radius - contact_distance - 0.01 * L() + max_coords = jnp.max(stent_vertices, axis=0) + max_stent_radius + influence_radius + contact_distance + 0.01 * L() mask = jnp.all((data_points >= min_coords) & (data_points <= max_coords), axis=1) return mask def compute_sdf_contact_displacements( data: dict, stent_vertices: jx.Array, s: float, - target_stent_radius: float, current_stent_radius: float, *, + target_stent_radius: float | jx.Array, current_stent_radius: float | jx.Array, *, influence_radius: float | None = None, contact_distance: float | None = None, - f_scale: float | None = None, + f_scale: float | None = None, cap_height_fraction: float = 1.0, ) -> tuple[np.ndarray, np.ndarray, float]: """Compute SDF-contact displacements for stent deployment. @@ -603,16 +662,28 @@ def compute_sdf_contact_displacements( Stent axis vertices, shape ``(V, 3)``. s : float Signed force scale. - target_stent_radius : float - Target stent radius for SDF computation. - current_stent_radius : float - Current deployment radius of the stent. + target_stent_radius : float or jx.Array + Target stent radius for SDF computation (scalar, or per-vertex radii + of shape ``(V,)``; only the maximum is used, for bounding-box + culling). + current_stent_radius : float or jx.Array + Current deployment radius of the stent: scalar, or per-vertex radii + of shape ``(V,)``, linearly interpolated along each capsule segment + so the stent can follow an arbitrary radius profile along the + centerline (e.g. a flared/trumpet end, a taper, or a local + narrowing). influence_radius : float Radial distance beyond the stent within which points are displaced. contact_distance : float Distance threshold for stent–wall contact. f_scale : float Force magnitude scaling factor. + cap_height_fraction : float + Axial semi-axis of the capsule end caps as a fraction of the local + radius; ``1.0`` (default) gives spherical caps, smaller values + (e.g. ``0.35``) flatten them into half ellipsoids, which is + recommended for flared or concave radius profiles (see + :func:`_tapered_capsule_segment_distances`). Returns ------- @@ -668,7 +739,7 @@ def compute_sdf_contact_displacements( # segments ensures C¹-continuous distance and direction fields. start_time = time.time() logger.debug(f"Total # surface and centerline points combined: {query_points.shape[0]}") - combined_final_dist_to_surface, combined_final_direction = smin_sdf_capsule_contact_sculpt(query_points, stent_vertices, current_stent_radius) + combined_final_dist_to_surface, combined_final_direction = smin_sdf_capsule_contact_sculpt(query_points, stent_vertices, current_stent_radius, cap_height_fraction) combined_final_dist_to_surface = np.array(combined_final_dist_to_surface) combined_final_direction = np.array(combined_final_direction) diff --git a/svmorph/core/geometry.py b/svmorph/core/geometry.py index 86c40782..a33cde48 100644 --- a/svmorph/core/geometry.py +++ b/svmorph/core/geometry.py @@ -125,3 +125,90 @@ def resample_stent_axis( new_vertices[:, dim] = np.interp(new_s, subsegment_s, subsegment_points[:, dim]) return jnp.array(new_vertices) + + +def _normalized_arc_positions(stent_vertices: np.ndarray) -> np.ndarray: + """Arc-length position of each stent axis vertex, normalized to [0, 1] from the first vertex.""" + points = np.asarray(stent_vertices, dtype=float) + if len(points) < 2: + return np.zeros(len(points)) + segment_lengths = np.linalg.norm(np.diff(points, axis=0), axis=1) + arc_positions = np.concatenate(([0.0], np.cumsum(segment_lengths))) + total_length = arc_positions[-1] + if total_length <= 0.0: + return np.zeros(len(points)) + return arc_positions / total_length + + +def stent_radius_profile(stent_vertices: np.ndarray, control_points) -> np.ndarray: + """Per-vertex stent radii from arbitrary (position, radius) control points. + + Builds a variable radius profile for the tapered capsule-chain stent SDF + (see :func:`svmorph.core.deformation.compute_sdf_contact_displacements`): + one radius per stent axis vertex, linearly interpolated between the + control points along the stent axis. + + Parameters + ---------- + stent_vertices : np.ndarray + Stent axis vertices, shape ``(V, 3)`` (e.g. from + :func:`resample_stent_axis`). + control_points : iterable of (float, float) + ``(position, radius)`` pairs, where *position* is the normalized + arc-length position in ``[0, 1]`` measured from the first axis + vertex. + + Returns + ------- + np.ndarray + Per-vertex stent radii, shape ``(V,)``. + """ + control_points = sorted((float(position), float(radius)) for position, radius in control_points) + if not control_points: + raise ValueError("At least one radius profile control point is required") + positions = [position for position, _ in control_points] + radii = [radius for _, radius in control_points] + return np.interp(_normalized_arc_positions(stent_vertices), positions, radii) + + +def flared_stent_radius_profile( + stent_vertices: np.ndarray, body_radius: float, flare_radius: float, + flare_length: float, flare_at_axis_start: bool = False, +) -> np.ndarray: + """Per-vertex stent radii for a stent with one flared (funnel/trumpet) end. + + The radius transitions from *body_radius* to *flare_radius* over + *flare_length* at one end of the stent axis with a smoothstep profile. + The result can be passed (scaled to the current deployment radius) as the + per-vertex ``current_stent_radius`` of + :func:`svmorph.core.deformation.compute_sdf_contact_displacements`. + + Parameters + ---------- + stent_vertices : np.ndarray + Stent axis vertices, shape ``(V, 3)``. + body_radius : float + Stent radius away from the flared end. + flare_radius : float + Stent radius at the tip of the flared end (may also be smaller than + *body_radius* for a tapered stent). + flare_length : float + Length of the radius transition, measured along the stent axis from + the flared end, in the same unit as the vertex coordinates. + flare_at_axis_start : bool + Flare the first-vertex end of the stent axis instead of the + last-vertex end. + + Returns + ------- + np.ndarray + Per-vertex stent radii, shape ``(V,)``. + """ + points = np.asarray(stent_vertices, dtype=float) + segment_lengths = np.linalg.norm(np.diff(points, axis=0), axis=1) + arc_positions = np.concatenate(([0.0], np.cumsum(segment_lengths))) + total_length = arc_positions[-1] + distance_from_flared_end = arc_positions if flare_at_axis_start else total_length - arc_positions + t = np.clip(1.0 - distance_from_flared_end / float(flare_length), 0.0, 1.0) + t = t * t * (3.0 - 2.0 * t) # smoothstep + return float(body_radius) + (float(flare_radius) - float(body_radius)) * t diff --git a/tests/test_deformation.py b/tests/test_deformation.py new file mode 100644 index 00000000..6e83d592 --- /dev/null +++ b/tests/test_deformation.py @@ -0,0 +1,92 @@ +"""Tests for svmorph.core.deformation – tapered capsule-chain SDF.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from svmorph.core.deformation import capsule_sdf, smin_sdf_capsule_contact_sculpt +from svmorph.core.units import set_unit_scale + + +@pytest.fixture(autouse=True) +def _unit_scale(): + set_unit_scale(1.0) + + +def _straight_axis(n: int = 11, length: float = 1.0) -> np.ndarray: + """Return n evenly-spaced stent axis vertices along the Z axis from 0 to length.""" + z = np.linspace(0.0, length, n) + return np.column_stack([np.zeros(n), np.zeros(n), z]) + + +def _surface_radius_along_ray(axis_point: np.ndarray, radii, cap_height_fraction: float = 1.0, + vertices: np.ndarray | None = None) -> float: + """Radial position of the SDF zero crossing on a +X ray from axis_point.""" + vertices = _straight_axis() if vertices is None else vertices + radial = np.linspace(0.01, 2.0, 400) + query = np.tile(np.asarray(axis_point, dtype=float), (len(radial), 1)) + query[:, 0] += radial + distances = np.asarray(capsule_sdf(query, vertices, radii, cap_height_fraction))[:, 0] + return float(radial[np.argmin(np.abs(distances))]) + + +def test_scalar_radius_matches_constant_profile(): + vertices = _straight_axis() + query = np.array([[0.35, 0.1, 0.4], [0.0, 0.0, -0.3], [0.5, -0.2, 1.2]]) + scalar = np.asarray(capsule_sdf(query, vertices, 0.25)) + per_vertex = np.asarray(capsule_sdf(query, vertices, np.full(len(vertices), 0.25))) + np.testing.assert_allclose(scalar, per_vertex, atol=1e-6) + + +def test_tapered_radius_is_interpolated_along_segment(): + # Single cone segment from radius 0.2 to 0.4: surface radius at the midpoint is 0.3 + vertices = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) + radii = np.array([0.2, 0.4]) + midpoint_radius = _surface_radius_along_ray([0.0, 0.0, 0.5], radii, vertices=vertices) + assert midpoint_radius == pytest.approx(0.3, abs=0.01) + + +def test_smin_sculpt_matches_capsule_sdf(): + vertices = _straight_axis() + radii = np.linspace(0.2, 0.5, len(vertices)) + query = np.array([[0.4, 0.0, 0.2], [0.1, 0.3, 0.9], [0.0, 0.0, 1.4]]) + sdf = np.asarray(capsule_sdf(query, vertices, radii, 0.5)) + sculpt_dist, sculpt_dir = smin_sdf_capsule_contact_sculpt(query, vertices, radii, 0.5) + np.testing.assert_allclose(sdf, np.asarray(sculpt_dist), atol=1e-6) + assert np.asarray(sculpt_dir).shape == (len(query), 3) + + +def test_default_cap_is_spherical(): + # A point on the axis beyond the end of a uniform capsule chain is on the surface + # at one radius past the end vertex (spherical end cap) + vertices = _straight_axis() + r = 0.3 + query = np.array([[0.0, 0.0, 1.0 + r]]) + distance = float(np.asarray(capsule_sdf(query, vertices, r))[0, 0]) + assert distance == pytest.approx(0.0, abs=1e-6) + + +def test_flattened_cap_height(): + # With a flattened cap the surface on the axis is cap_height_fraction * radius + # past the end vertex + vertices = _straight_axis() + r, fraction = 0.3, 0.35 + query = np.array([[0.0, 0.0, 1.0 + fraction * r]]) + distance = float(np.asarray(capsule_sdf(query, vertices, r, fraction))[0, 0]) + assert distance == pytest.approx(0.0, abs=1e-6) + + +def test_concave_profile_preserved_by_flattened_caps(): + # Dumbbell profile: wide - narrow waist - wide. With spherical caps the wide + # capsules' end caps bulge into the waist and fill it in; flattened caps keep it. + vertices = _straight_axis(n=21, length=2.0) + radii = np.full(len(vertices), 1.2) + radii[6:15] = 0.5 # waist between z=0.6 and z=1.4 + waist_center = [0.0, 0.0, 1.0] + flattened_waist = _surface_radius_along_ray(waist_center, radii, cap_height_fraction=0.35, + vertices=vertices) + spherical_waist = _surface_radius_along_ray(waist_center, radii, cap_height_fraction=1.0, + vertices=vertices) + assert flattened_waist == pytest.approx(0.5, abs=0.15) + assert spherical_waist > 1.0 diff --git a/tests/test_geometry.py b/tests/test_geometry.py index f7bf14ce..1e1c9bf9 100644 --- a/tests/test_geometry.py +++ b/tests/test_geometry.py @@ -1,9 +1,13 @@ -"""Tests for svmorph.core.geometry – centerline resampling.""" +"""Tests for svmorph.core.geometry – centerline resampling and radius profiles.""" import numpy as np import pytest -from svmorph.core.geometry import resample_stent_axis +from svmorph.core.geometry import ( + flared_stent_radius_profile, + resample_stent_axis, + stent_radius_profile, +) def _straight_line(n: int = 20, length: float = 2.0) -> np.ndarray: @@ -70,3 +74,39 @@ def test_too_few_points_raises(): ptm, sbm = _flat_maps(1) with pytest.raises(ValueError, match="Not enough points"): resample_stent_axis(pts, ptm, sbm, 0, 0.5, 0.1, sampling_direction=-1) + + +def test_stent_radius_profile_interpolates_control_points(): + pts = _straight_line(n=11, length=1.0) + radii = stent_radius_profile(pts, [(0.0, 0.2), (0.5, 0.4), (1.0, 0.2)]) + assert radii.shape == (11,) + assert radii[0] == pytest.approx(0.2) + assert radii[5] == pytest.approx(0.4) + assert radii[-1] == pytest.approx(0.2) + assert radii[2] == pytest.approx(0.2 + 0.2 * (0.2 / 0.5), abs=1e-6) + + +def test_stent_radius_profile_requires_control_points(): + pts = _straight_line(n=5) + with pytest.raises(ValueError, match="control point"): + stent_radius_profile(pts, []) + + +def test_flared_stent_radius_profile_at_axis_end(): + pts = _straight_line(n=21, length=2.0) + radii = flared_stent_radius_profile(pts, body_radius=0.4, flare_radius=0.6, flare_length=0.5) + assert radii.shape == (21,) + # Body radius away from the flare, flare radius exactly at the flared (last) vertex + assert radii[0] == pytest.approx(0.4) + assert radii[10] == pytest.approx(0.4) + assert radii[-1] == pytest.approx(0.6) + # Monotonic smoothstep transition within the flare length + transition = radii[15:] + assert np.all(np.diff(transition) >= 0.0) + + +def test_flared_stent_radius_profile_at_axis_start(): + pts = _straight_line(n=21, length=2.0) + radii = flared_stent_radius_profile(pts, 0.4, 0.6, 0.5, flare_at_axis_start=True) + assert radii[0] == pytest.approx(0.6) + assert radii[-1] == pytest.approx(0.4) From 7aaaa5d2441b5dfcdf8edbb3ff4d5c83544e9925 Mon Sep 17 00:00:00 2001 From: Andras Lasso Date: Wed, 9 Sep 2026 23:19:36 -0400 Subject: [PATCH 2/2] ENH: Add flared stent options to the deploy_stent script New command line options for the variable-radius SDFStent deployment: --flare-R, --flare-length and --flare-end deploy a stent with a flared (funnel/trumpet shaped) end whose radius transitions from the target radius to the flare radius with a smoothstep profile; --cap-height-fraction flattens the capsule end caps into half ellipsoids (recommended ~0.35 for flared or concave profiles; the default 1.0 keeps the classic spherical caps). During deployment the whole radius profile is expanded proportionally with the nominal radius, so the flared end reaches the flare radius exactly when the body reaches the target radius. Documented in the README with a usage example. Co-Authored-By: Claude Fable 5 --- README.md | 22 ++++++++++++++ svmorph/scripts/deploy_stent.py | 51 +++++++++++++++++++++++++++++++-- 2 files changed, 71 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index ddbb9550..5eda1b42 100644 --- a/README.md +++ b/README.md @@ -295,6 +295,28 @@ python -m svmorph.scripts.deploy_stent \ --out-mesh deployed_surface.vtp --out-cl deployed_centerline.vtp ``` +**Deploy a flared stent (SDFStent, variable radius)** — same as above, but the +stent radius follows a profile along the centerline instead of being uniform: +the end opposite the `--start` point flares from the 0.4 cm body radius to +0.6 cm over the last 0.5 cm of the stent (funnel/trumpet shape). +`--cap-height-fraction` flattens the capsule end caps into half ellipsoids +(axial semi-axis = fraction × local radius); this keeps the wide flared end +from inflating a spherical bulge into the vessel beyond the stent end, and +also lets concave radius profiles (local narrowings) deploy without being +filled in by the neighboring wide capsules. Arbitrary radius profiles can be +built programmatically with `svmorph.core.stent_radius_profile` (control +points) or `svmorph.core.flared_stent_radius_profile` and passed as the +per-vertex `current_stent_radius` of `compute_sdf_contact_displacements`. + +```bash +python -m svmorph.scripts.deploy_stent \ + --mesh surface.vtp --cline centerline.vtp \ + --start 123 --target-R 0.4 --start-R 0.05 --length 1.7 \ + --flare-R 0.6 --flare-length 0.5 --flare-end end \ + --cap-height-fraction 0.35 \ + --out-mesh deployed_surface.vtp --out-cl deployed_centerline.vtp +``` + **Deploy with concurrent axis straightening (SDFStent)** — same stent geometry as above, but after each expansion step the stent axis is projected toward the straight line connecting its endpoints (strength 0.075), gradually removing curvature diff --git a/svmorph/scripts/deploy_stent.py b/svmorph/scripts/deploy_stent.py index 13ef848c..94bcf80e 100644 --- a/svmorph/scripts/deploy_stent.py +++ b/svmorph/scripts/deploy_stent.py @@ -24,6 +24,8 @@ import argparse import time +import numpy as np + from svmorph.core import deformation, geometry, mesh_data from svmorph.core.units import L from svmorph.logging import get_logger @@ -45,6 +47,26 @@ def build_parser() -> argparse.ArgumentParser: parser.add_argument("--target-R", type=float, default=None, help="Target deployed stent radius (default: 0.4 cm)") parser.add_argument("--start-R", type=float, default=None, help="Initial crimped stent radius (default: 0.05 cm)") parser.add_argument("--length", type=float, default=None, help="Stent length along centerline (default: 1.7 cm)") + parser.add_argument( + "--flare-R", type=float, default=None, + help="Optional flared (funnel/trumpet) end: stent radius at the tip of the flared end " + "(default: no flare)", + ) + parser.add_argument( + "--flare-length", type=float, default=None, + help="Length of the flare radius transition along the stent axis (default: 0.5 cm)", + ) + parser.add_argument( + "--flare-end", choices=("start", "end"), default="end", + help="Which end of the stent axis is flared: 'start' is the --start point side, " + "'end' is the opposite (proximally resampled) side", + ) + parser.add_argument( + "--cap-height-fraction", type=float, default=1.0, + help="Axial height of the capsule end caps as a fraction of the local stent radius; " + "1.0 keeps the classic spherical caps, smaller values (e.g. 0.35) flatten them " + "into half ellipsoids, recommended for flared or concave radius profiles", + ) parser.set_defaults(out_mesh="deployed_surface.vtp", out_cl="deployed_centerline.vtp") return parser @@ -80,6 +102,22 @@ def main(argv: list[str] | None = None) -> None: ) logger.info(f"Stent axis: {len(axis_pts)} vertices over {deployed_length:.2f}") + # Optional flared end: per-vertex radius profile as fractions of the nominal target + # radius, so the whole profile can be scaled proportionally during deployment + radius_profile_fractions = None + if args.flare_R is not None: + if args.flare_length is None: + args.flare_length = 0.5 * L() + target_radii = geometry.flared_stent_radius_profile( + np.asarray(axis_pts), args.target_R, args.flare_R, args.flare_length, + flare_at_axis_start=(args.flare_end == "start"), + ) + radius_profile_fractions = target_radii / args.target_R + logger.info( + f"Flared '{args.flare_end}' end: radius {args.target_R:.4f} -> {args.flare_R:.4f} " + f"over {args.flare_length:.2f}" + ) + a, b = mesh_data.compute_material_constants(1.0, 0.2) snapshots = common.SnapshotManager( @@ -99,12 +137,21 @@ def main(argv: list[str] | None = None) -> None: iteration = 0 t0 = time.time() while True: + if radius_profile_fractions is None: + current_stent_radius = cur_R + bounding_stent_radius = args.target_R + else: + # Scale the whole radius profile proportionally with the nominal radius, keeping + # the smooth-min smoothing offset constant along the stent + current_stent_radius = radius_profile_fractions * (cur_R + smoothing_k) - smoothing_k + bounding_stent_radius = float(radius_profile_fractions.max()) * args.target_R surf_disp, cl_disp, dR = deformation.compute_sdf_contact_displacements( ctx.data, axis_pts, s=-1.0, - target_stent_radius=args.target_R, - current_stent_radius=cur_R, + target_stent_radius=bounding_stent_radius, + current_stent_radius=current_stent_radius, + cap_height_fraction=args.cap_height_fraction, ) if cur_R + dR > args.target_R: logger.info("Next increment would overshoot target -- done.")