Upload endograft_script.py

endograft_script.py

@@ -2,104 +2,107 @@ import vtk
 import numpy as np
 import tempfile
 
-def arc_length(centerline):
-    arc_lengths = np.zeros(len(centerline))
-    for i in range(1, len(centerline)):
-        arc_lengths[i] = arc_lengths[i-1] + np.linalg.norm(centerline[i] - centerline[i-1])
-    return arc_lengths
+def arc_length(cl):
+    s = np.zeros(len(cl))
+    for i in range(1, len(cl)):
+        s[i] = s[i-1] + np.linalg.norm(cl[i] - cl[i-1])
+    return s
 
 def compute_parallel_frames(centerline):
     frames = []
 
-    tangent_0 = centerline[1] - centerline[0]
-    tangent_0 /= np.linalg.norm(tangent_0)
+    t0 = centerline[1] - centerline[0]
+    t0 /= np.linalg.norm(t0)
 
-    reference_vector = np.array([0, 0, 1])     # selecting z-axis as reference for calculating initial normal 
-    if abs(np.dot(tangent_0, reference_vector)) > 0.9:
-        reference_vector = np.array([1, 0, 0])
+    ref = np.array([0, 0, 1])
+    if abs(np.dot(t0, ref)) > 0.9:
+        ref = np.array([1, 0, 0])
 
-    normal_0 = np.cross(tangent_0, reference_vector)
-    normal_0 /= np.linalg.norm(normal_0)
-    binormal_0 = np.cross(tangent_0, normal_0)
+    n0 = np.cross(t0, ref)
+    n0 /= np.linalg.norm(n0)
+    b0 = np.cross(t0, n0)
 
-    frames.append((tangent_0, normal_0, binormal_0))
+    frames.append((t0, n0, b0))
 
-    for point_idx in range(1, len(centerline)):
-        tangent_point = centerline[point_idx] - centerline[point_idx-1]
-        tangent_point /= np.linalg.norm(tangent_point)
+    for i in range(1, len(centerline)):
+        ti = centerline[i] - centerline[i-1]
+        ti /= np.linalg.norm(ti)
 
-        tangent_prev, normal_prev, binormal_prev = frames[-1]
-        rotation_axis = np.cross(tangent_prev, tangent_point)
+        t_prev, n_prev, b_prev = frames[-1]
+        v = np.cross(t_prev, ti)
 
-        if np.linalg.norm(rotation_axis) < 1e-6:
-            frames.append((tangent_point, normal_prev, binormal_prev))
+        if np.linalg.norm(v) < 1e-6:
+            frames.append((ti, n_prev, b_prev))
             continue
 
-        rotation_axis /= np.linalg.norm(rotation_axis)
-        angle = np.arccos(np.clip(np.dot(tangent_prev, tangent_point), -1, 1))
+        v /= np.linalg.norm(v)
+        angle = np.arccos(np.clip(np.dot(t_prev, ti), -1, 1))
 
-        def rotate_vector(vector):
-            return (vector * np.cos(angle) + np.cross(rotation_axis, vector) * np.sin(angle) + rotation_axis * np.dot(rotation_axis, vector) * (1 - np.cos(angle)))
+        def rot(vec):
+            return (vec * np.cos(angle) + np.cross(v, vec) * np.sin(angle) + v * np.dot(v, vec) * (1 - np.cos(angle)))
 
-        frames.append((tangent_point, rotate_vector(normal_prev), rotate_vector(binormal_prev)))
+        frames.append((ti, rot(n_prev), rot(b_prev)))
 
     return frames
 
-def project_point_to_surface(point, locator):
-    
+def project_point_to_surface(p, surface_poly):
+    locator = vtk.vtkCellLocator()
+    locator.SetDataSet(surface_poly)
+    locator.BuildLocator()
+
     closest = [0.0, 0.0, 0.0]
     cellId = vtk.mutable(0)
     subId = vtk.mutable(0)
     dist2 = vtk.mutable(0.0)
 
-    locator.FindClosestPoint(point, closest, cellId, subId, dist2)
+    locator.FindClosestPoint(p, closest, cellId, subId, dist2)
     return np.array(closest)
 
-def sample_centerline_by_distance(centerline, arc_lengths, frames, query_distance):
+def sample_centerline_by_s(centerline, s_cl, frames, s):
+    idx = np.searchsorted(s_cl, s)
 
-    idx = np.searchsorted(arc_lengths, query_distance)
     if idx <= 0:
         return centerline[0], frames[0]
     if idx >= len(centerline):
         return centerline[-1], frames[-1]
 
-    arc_length_0, arc_length_1 = arc_lengths[idx-1], arc_lengths[idx]
-    point_0, point_1 = centerline[idx-1], centerline[idx]
-    tangent = (query_distance - arc_length_0) / (arc_length_1 - arc_length_0 + 1e-12)
+    s0, s1 = s_cl[idx-1], s_cl[idx]
+    P0, P1 = centerline[idx-1], centerline[idx]
+    t = (s - s0) / (s1 - s0 + 1e-12)
 
-    point = (1 - tangent) * point_0 + tangent * point_1
+    P = (1 - t) * P0 + t * P1
 
-    tangent_0, normal_0, binormal_0 = frames[idx-1]
-    tangent_1, normal_1, binormal_1 = frames[idx]
+    t0, N0, B0 = frames[idx-1]
+    t1, N1, B1 = frames[idx]
 
-    tangent_vector = (1 - tangent) * tangent_0 + tangent * tangent_1
-    normal_vector = (1 - tangent) * normal_0 + tangent * normal_1
-    binormal_vector = (1 - tangent) * binormal_0 + tangent * binormal_1
+    t_vec = (1 - t) * t0 + t * t1
+    N_vec = (1 - t) * N0 + t * N1
+    B_vec = (1 - t) * B0 + t * B1
 
-    tangent_vector /= np.linalg.norm(tangent_vector)
-    normal_vector /= np.linalg.norm(normal_vector)
-    binormal_vector /= np.linalg.norm(binormal_vector)
+    t_vec /= np.linalg.norm(t_vec)
+    N_vec /= np.linalg.norm(N_vec)
+    B_vec /= np.linalg.norm(B_vec)
 
-    return point, (tangent_vector, normal_vector, binormal_vector)
+    return P, (t_vec, N_vec, B_vec)
 
-# Fabric Component
-def create_component_fabric(centerline, frames, start_position, length, diameter, start_amplitude_pp=2.0, waves=8):
+# Fabric
+def create_component_fabric(centerline, frames, start_s, length, diameter, start_amp_pp=2.0, waves=8):
     # shift centerline sampling
     shifted_centerline = centerline.copy()
 
-    arc_lengths = arc_length(centerline)
+    s_cl = arc_length(centerline)
 
     pts = vtk.vtkPoints()
     polys = vtk.vtkCellArray()
     rings = []
 
-    A = start_amplitude_pp / 2.0
+    A = start_amp_pp / 2.0
     n_sections = 180
     n_theta = 64
 
     for i in range(n_sections):
-        distance_base = start_position + length * i / (n_sections - 1)
-        point_base, (tangent, normal, binormal) = sample_centerline_by_distance(centerline, arc_lengths, frames, distance_base)
+        s_base = start_s + length * i / (n_sections - 1)
+        P_base, (t, N, B) = sample_centerline_by_s(centerline, s_cl, frames, s_base)
 
         ring = []
 
@@ -107,14 +110,14 @@ def create_component_fabric(centerline, frames, start_position, length, diameter
             theta = 2 * np.pi * j / n_theta
             axial_offset = 0.0
 
-            if i == 0 and start_amplitude_pp > 0:
+            if i == 0 and start_amp_pp > 0:
                 axial_offset = A * np.sin(waves * theta) - A
 
-            point = point_base + axial_offset * tangent
+            P = P_base + axial_offset * t
             r = diameter / 2.0
-            position = point + r * (np.cos(theta) * normal + np.sin(theta) * binormal)
+            pos = P + r * (np.cos(theta) * N + np.sin(theta) * B)
 
-            ring.append(pts.InsertNextPoint(position))
+            ring.append(pts.InsertNextPoint(pos))
 
         rings.append(ring)
 
@@ -135,13 +138,14 @@ def create_component_fabric(centerline, frames, start_position, length, diameter
     
 # Struts
 
-def create_proximal_bare_stent(centerline, frames, fabric_start_distance, main_diam_mm, length, main_fabric, peak_to_peak_amplitude=15.0, 
-                               peak_spacing=7.0, component_type="proximal", wire_radius=0.35, waves=8, samples=240):
+def create_proximal_bare_stent(centerline, frames, fabric_start_s, main_diam_mm, length, main_fabric, peak_to_peak=15.0, 
+                               peak_spacing=7.0, component_type="proximal", wire_radius=0.35, waves=8, samples=240
+                               ):
 
-    arc_lengths = arc_length(centerline)
+    s_cl = arc_length(centerline)
     append = vtk.vtkAppendPolyData()
 
-    # build fabric normals & locator once
+    # --- build fabric normals & locator ONCE ---
     normal_gen = vtk.vtkPolyDataNormals()
     normal_gen.SetInputData(main_fabric)
     normal_gen.ComputePointNormalsOn()
@@ -150,35 +154,38 @@ def create_proximal_bare_stent(centerline, frames, fabric_start_distance, main_d
     normal_gen.Update()
     fabric_normals = normal_gen.GetOutput().GetPointData().GetNormals()
 
-    point_locator = vtk.vtkPointLocator()
-    point_locator.SetDataSet(main_fabric)
-    point_locator.BuildLocator()
+    locator = vtk.vtkPointLocator()
+    locator.SetDataSet(main_fabric)
+    locator.BuildLocator()
+    # -----------------------------------------
 
-    A = peak_to_peak_amplitude              # peak-peak amplitude = 15
+    A = peak_to_peak              # peak-peak amplitude = 15
     r = main_diam_mm / 2.0
     gap = peak_spacing
     step = A + gap
 
     n_crowns = int(np.ceil(length / step))
-    distance_peaks = [fabric_start_distance + 1.0 + i * step for i in range(n_crowns)]
-    if ((component_type == "distal" or component_type == "extension") and (len(distance_peaks) > 4)):
-      distance_peaks = distance_peaks[3:]
-
-# fication strut 
-    A_extra = peak_to_peak_amplitude / 2.0      # 7.5 mm height
+    s_peaks = [fabric_start_s + 1.0 + i * step for i in range(n_crowns)]
+    if ((component_type == "distal" or component_type == "extension") and (len(s_peaks) > 4)):
+      s_peaks = s_peaks[3:]
+
+# =========================================================
+# ===== EXTRA CLIPPED HALF sin =================
+# =========================================================
+    A_extra = peak_to_peak / 2.0      # 7.5 mm height
     extra_top = -1e-10                  # peak at -7.5
 
     if component_type == "proximal":
-      point_0 = centerline[0]
-      tangent_0, normal_0, binormal_0 = frames[0]
-      point = point_0 + extra_top * tangent_0   # extrapolate backward
-      tangent, normal, binormal = tangent_0, normal_0, binormal_0
+      P0 = centerline[0]
+      t0, N0, B0 = frames[0]
+      P = P0 + extra_top * t0   # extrapolate backward
+      t, N, B = t0, N0, B0
 
     elif component_type == "distal":
-      point, (tangent, normal, binormal) = sample_centerline_by_distance(centerline, arc_lengths, frames, fabric_start_distance + length - extra_top)
+      P, (t, N, B) = sample_centerline_by_s(centerline, s_cl, frames, fabric_start_s + length - extra_top)
 
     else:
-        point, (tangent, normal, binormal) = sample_centerline_by_distance(centerline, arc_lengths, frames, extra_top)
+        P, (t, N, B) = sample_centerline_by_s(centerline, s_cl, frames, extra_top)
 
     pts = vtk.vtkPoints()
     lines = vtk.vtkCellArray()
@@ -189,22 +196,22 @@ def create_proximal_bare_stent(centerline, frames, fabric_start_distance, main_d
         # keep ONLY negative or positive half of sine
         if component_type == "proximal":
               axial_offset = A_extra * min(0.0, np.sin(waves * theta))
-              raw_position = (point
-                  + (r + wire_radius) * (np.cos(theta)*normal + np.sin(theta)*binormal)
-                  + axial_offset * tangent)
+              raw_pos = (P
+                  + (r + wire_radius) * (np.cos(theta)*N + np.sin(theta)*B)
+                  + axial_offset * t)
         elif (component_type == "distal" and length != 0):
               axial_offset = A_extra * max(0.0, np.sin(waves * theta))
 
-              raw_position = (point
-                        + (r + wire_radius) * (np.cos(theta)*normal + np.sin(theta)*binormal)
-                        + axial_offset * tangent)
+              raw_pos = (P
+                        + (r + wire_radius) * (np.cos(theta)*N + np.sin(theta)*B)
+                        + axial_offset * t)
         else:
           continue
 
-        radial = (np.cos(theta)*normal + np.sin(theta)*binormal)
+        radial = (np.cos(theta)*N + np.sin(theta)*B)
         radial /= np.linalg.norm(radial)
 
-        pos = point + (r + wire_radius) * radial + axial_offset * tangent
+        pos = P + (r + wire_radius) * radial + axial_offset * t
         pts.InsertNextPoint(pos)
 
         if i > 0:
@@ -225,17 +232,14 @@ def create_proximal_bare_stent(centerline, frames, fabric_start_distance, main_d
     tube.Update()
 
     append.AddInputData(tube.GetOutput())
-    
-    # normal struts along the graft body
-
-    cell_locator = vtk.vtkCellLocator()
-    cell_locator.SetDataSet(main_fabric)
-    cell_locator.BuildLocator()
+    # =========================================================
+    # ================= NORMAL FULL struts ====================
+    # =========================================================
 
-    for distance_peak in distance_peaks:
-        point, (tangent, normal, binormal) = sample_centerline_by_distance(centerline, arc_lengths, frames, distance_peak)
+    for s0 in s_peaks:
+        P, (t, N, B) = sample_centerline_by_s(centerline, s_cl, frames, s0)
 
-        points = vtk.vtkPoints()
+        pts = vtk.vtkPoints()
         lines = vtk.vtkCellArray()
 
         for i in range(samples + 1):
@@ -243,19 +247,19 @@ def create_proximal_bare_stent(centerline, frames, fabric_start_distance, main_d
 
             axial_offset = A * (1 + np.sin(waves * theta)) / 2.0
 
-            raw_position = (point
-                       + (r + wire_radius) * (np.cos(theta)*normal + np.sin(theta)*binormal)
-                       + axial_offset * tangent)
+            raw_pos = (P
+                       + (r + wire_radius) * (np.cos(theta)*N + np.sin(theta)*B)
+                       + axial_offset * t)
 
-            surface_point = project_point_to_surface(raw_position, cell_locator)
+            surf_p = project_point_to_surface(raw_pos, main_fabric)
 
-            pid = point_locator.FindClosestPoint(surface_point)
-            normal = np.array(fabric_normals.GetTuple(pid))
-            normal /= np.linalg.norm(normal)
+            pid = locator.FindClosestPoint(surf_p)
+            n = np.array(fabric_normals.GetTuple(pid))
+            n /= np.linalg.norm(n)
 
-            position = surface_point + wire_radius * normal
+            pos = surf_p + wire_radius * n
 
-            points.InsertNextPoint(position)
+            pts.InsertNextPoint(pos)
 
             if i > 0:
                 l = vtk.vtkLine()
@@ -264,7 +268,7 @@ def create_proximal_bare_stent(centerline, frames, fabric_start_distance, main_d
                 lines.InsertNextCell(l)
 
         poly = vtk.vtkPolyData()
-        poly.SetPoints(points)
+        poly.SetPoints(pts)
         poly.SetLines(lines)
 
         tube = vtk.vtkTubeFilter()
@@ -279,9 +283,9 @@ def create_proximal_bare_stent(centerline, frames, fabric_start_distance, main_d
     append.Update()
     return append.GetOutput()
 
-def create_component_stents(centerline, frames, fabric_start_distance, length, diameter, fabric, peak_to_peak=15.0, peak_spacing=7.0, component_type="proximal"):
+def create_component_stents(centerline, frames, fabric_start_s, length, diameter, fabric, peak_to_peak=15.0, peak_spacing=7.0, component_type="proximal"):
 
-    stents = create_proximal_bare_stent(centerline, frames, fabric_start_distance, diameter, length, fabric, peak_to_peak, peak_spacing, component_type)
+    stents = create_proximal_bare_stent(centerline, frames, fabric_start_s, diameter, length, fabric, peak_to_peak, peak_spacing, component_type)
 
     return stents
 
@@ -299,7 +303,7 @@ def merge_components(proximal, extension, distal):
 
     return clean
 
-def create_endograft(centerline, diameter, proximal_length, extension_length, distal_length, proximal_gap, extension_gap, distal_gap):
+def create_endograft(centerline, diameter, prox_length,extension_length, distal_length, prox_gap, extension_gap, distal_gap):
     frames = compute_parallel_frames(centerline)
     peak_to_peak = 15.0
 
@@ -307,11 +311,11 @@ def create_endograft(centerline, diameter, proximal_length, extension_length, di
     # PROXIMAL COMPONENT
     # -------------------------
 
-    proximal_start = 0.0
+    prox_start = 0.0
 
-    proximal_fabric = create_component_fabric(centerline, frames, proximal_start, proximal_length, diameter, start_amplitude_pp=2.0)
+    prox_fabric = create_component_fabric(centerline, frames, prox_start, prox_length, diameter, start_amp_pp=2.0)
 
-    proximal_stents = create_component_stents(centerline, frames, proximal_start, proximal_length, diameter, proximal_fabric, peak_spacing=proximal_gap, component_type="proximal")
+    prox_stents = create_component_stents(centerline, frames, prox_start, prox_length, diameter, prox_fabric, peak_spacing=prox_gap, component_type="proximal")
 
 
     # -------------------------
@@ -322,8 +326,8 @@ def create_endograft(centerline, diameter, proximal_length, extension_length, di
     else:
         prox_extension_overlap = 2*peak_to_peak + extension_gap
 
-    extension_start = proximal_length - prox_extension_overlap
-    extension_fabric = create_component_fabric(centerline, frames, extension_start, extension_length, diameter, start_amplitude_pp=0.0)
+    extension_start = prox_length - prox_extension_overlap
+    extension_fabric = create_component_fabric(centerline, frames, extension_start, extension_length, diameter, start_amp_pp=0.0)
 
     extension_stents = create_component_stents(centerline, frames, extension_start, extension_length, diameter, extension_fabric, peak_spacing=extension_gap, component_type="extension")
 
@@ -331,14 +335,14 @@ def create_endograft(centerline, diameter, proximal_length, extension_length, di
     # DISTAL COMPONENT
     # -------------------------
     distal_overlap = 3*peak_to_peak + 2*distal_gap
-    distal_start = proximal_length + extension_length - prox_extension_overlap - distal_overlap
-    distal_fabric = create_component_fabric(centerline, frames, distal_start, distal_length, diameter, start_amplitude_pp=0.0)
+    distal_start = prox_length + extension_length - prox_extension_overlap - distal_overlap
+    distal_fabric = create_component_fabric(centerline, frames, distal_start, distal_length, diameter, start_amp_pp=0.0)
     distal_stents = create_component_stents(centerline, frames, distal_start, distal_length, diameter, distal_fabric, peak_spacing=distal_gap, component_type="distal")
     # -------------------------
     # MERGE COMPONENTS
     # -------------------------
-    fabric_part = merge_components(proximal_fabric, extension_fabric, distal_fabric)
-    stent_part = merge_components(proximal_stents, extension_stents, distal_stents)
+    fabric_part = merge_components(prox_fabric, extension_fabric, distal_fabric)
+    stent_part = merge_components(prox_stents, extension_stents, distal_stents)
 
     endograft = vtk.vtkAppendPolyData()