Upload endograft_generation.py

endograft_generation.py

@@ -9,218 +9,98 @@ def arc_length(centerline):
         arc_lengths[i] = arc_lengths[i-1] + np.linalg.norm(centerline[i] - centerline[i-1])
     return arc_lengths
 
-def smooth_centerline(centerline, smoothing=0.5, samples=500):
-    centerline = np.asarray(centerline)
-    
-    # 1. Fit the initial spline
-    spline, _ = splprep(centerline.T, s=smoothing, k=3)
-    
-    # 2. Evaluate at extremely high resolution to measure true physical arc-length
-    u_fine = np.linspace(0, 1, 3000)
-    fine_points = np.vstack(splev(u_fine, spline)).T
-    
-    fine_arc_lengths = np.zeros(len(fine_points))
-    for i in range(1, len(fine_points)):
-        fine_arc_lengths[i] = fine_arc_lengths[i-1] + np.linalg.norm(fine_points[i] - fine_points[i-1])
-        
-    # 3. Re-sample the points so they are STRICTLY equidistant in true 3D space
-    target_arc_lengths = np.linspace(0, fine_arc_lengths[-1], samples)
-    
-    smoothed_centerline = np.zeros((samples, 3))
-    for i in range(3):
-        smoothed_centerline[:, i] = np.interp(target_arc_lengths, fine_arc_lengths, fine_points[:, i])
-        
-    return smoothed_centerline
-
 def compute_parallel_frames(centerline):
-    """
-    Uses Wang's Double Reflection Rotation Minimizing Frame (RMF) algorithm.
-    Mathematically guarantees zero twist along the 3D curve, eliminating strut shear.
-    """
     frames = []
-    
-    # Initial frame
-    t0 = centerline[1] - centerline[0]
-    t0 /= np.linalg.norm(t0)
-    
-    ref = np.array([0.0, 0.0, 1.0])
-    if abs(np.dot(t0, ref)) > 0.9:
-        ref = np.array([1.0, 0.0, 0.0])
-        
-    n0 = np.cross(t0, ref)
-    n0 /= np.linalg.norm(n0)
-    b0 = np.cross(t0, n0)
-    
-    frames.append((t0, n0, b0))
-    
-    for i in range(1, len(centerline) - 1):
-        x0 = centerline[i-1]
-        x1 = centerline[i]
-        x2 = centerline[i+1]
-        
-        t1 = x2 - x0
-        t1 /= np.linalg.norm(t1)
-        
-        t0, n0, b0 = frames[-1]
-        
-        # Wang's Double Reflection Method
-        v1 = x1 - x0
-        c1 = np.dot(v1, v1)
-        if c1 == 0:
-            frames.append((t1, n0, b0))
-            continue
-            
-        rL_t0 = t0 - (2.0 / c1) * np.dot(v1, t0) * v1
-        rL_n0 = n0 - (2.0 / c1) * np.dot(v1, n0) * v1
-        
-        v2 = t1 - rL_t0
-        c2 = np.dot(v2, v2)
-        if c2 == 0:
-            n1 = rL_n0
-        else:
-            n1 = rL_n0 - (2.0 / c2) * np.dot(v2, rL_n0) * v2
-            
-        n1 /= np.linalg.norm(n1)
-        b1 = np.cross(t1, n1)
-        
-        frames.append((t1, n1, b1))
-        
-    # Append the last frame 
-    t_last = centerline[-1] - centerline[-2]
-    t_last /= np.linalg.norm(t_last)
-    _, n_last, b_last = frames[-1]
-    frames.append((t_last, n_last, b_last))
-    
-    return frames
 
-def sample_centerline_by_distance(centerline, arc_lengths, frames, query_distance):
-    # Safe linear extrapolation for fixation struts that extend past the centerline
-    if query_distance <= 0:
-        t0, n0, b0 = frames[0]
-        return centerline[0] + query_distance * t0, (t0, n0, b0)
-    if query_distance >= arc_lengths[-1]:
-        t1, n1, b1 = frames[-1]
-        return centerline[-1] + (query_distance - arc_lengths[-1]) * t1, (t1, n1, b1)
+    tangent_0 = centerline[1] - centerline[0]
+    tangent_0 /= np.linalg.norm(tangent_0)
 
-    idx = np.searchsorted(arc_lengths, query_distance)
-    arc_length_0, arc_length_1 = arc_lengths[idx-1], arc_lengths[idx]
-    tangent_weight = (query_distance - arc_length_0) / (arc_length_1 - arc_length_0 + 1e-12)
+    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])
 
-    point = (1 - tangent_weight) * centerline[idx-1] + tangent_weight * centerline[idx]
+    normal_0 = np.cross(tangent_0, reference_vector)
+    normal_0 /= np.linalg.norm(normal_0)
+    binormal_0 = np.cross(tangent_0, normal_0)
 
-    t0, n0, b0 = frames[idx-1]
-    t1, n1, b1 = frames[idx]
+    frames.append((tangent_0, normal_0, binormal_0))
 
-    t_vec = (1 - tangent_weight) * t0 + tangent_weight * t1
-    n_vec = (1 - tangent_weight) * n0 + tangent_weight * n1
-    b_vec = (1 - tangent_weight) * b0 + tangent_weight * b1
+    for point_idx in range(1, len(centerline)):
+        tangent_point = centerline[point_idx] - centerline[point_idx-1]
+        tangent_point /= np.linalg.norm(tangent_point)
 
-    t_vec /= np.linalg.norm(t_vec)
-    n_vec /= np.linalg.norm(n_vec)
-    b_vec /= np.linalg.norm(b_vec)
+        tangent_prev, normal_prev, binormal_prev = frames[-1]
+        rotation_axis = np.cross(tangent_prev, tangent_point)
 
-    return point, (t_vec, n_vec, b_vec)
+        if np.linalg.norm(rotation_axis) < 1e-6:
+            frames.append((tangent_point, normal_prev, binormal_prev))
+            continue
 
+        rotation_axis /= np.linalg.norm(rotation_axis)
+        angle = np.arccos(np.clip(np.dot(tangent_prev, tangent_point), -1, 1))
 
-def create_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, phase_angle=0.0):
+        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)))
 
-    arc_lengths = arc_length(centerline)
-    phase_rad = np.radians(phase_angle)
-    append = vtk.vtkAppendPolyData()
+        frames.append((tangent_point, rotate_vector(normal_prev), rotate_vector(binormal_prev)))
 
-    A = peak_to_peak_amplitude
-    r = main_diam_mm / 2.0
-    gap = peak_spacing
-    step = A + gap
+    return frames
 
-    n_crowns = int(np.ceil(length / step))
-    distance_peaks = [fabric_start_distance + 1.0 + i * step for i in range(n_crowns)]
+def project_point_to_surface(point, locator, polydata):
     
-    if ((component_type == "distal" or component_type == "extension") and (len(distance_peaks) > 4)):
-        distance_peaks = distance_peaks[3:]
-    elif (component_type =="extension" and (2 < len(distance_peaks) <= 4)):
-        distance_peaks = distance_peaks[2:]
-
-    # ---------------------------------------------------------
-    # 1. Fixation Struts
-    # ---------------------------------------------------------
-    A_extra = peak_to_peak_amplitude / 2.0
-    pts = vtk.vtkPoints()
-    lines = vtk.vtkCellArray()
+    closest = [0.0, 0.0, 0.0]
+    cellId = vtk.mutable(0)
+    subId = vtk.mutable(0)
+    dist2 = vtk.mutable(0.0)
 
-    for i in range(samples + 1):
-        base_theta = 2 * np.pi * i / samples
-        theta_rot = base_theta + phase_rad
+    locator.FindClosestPoint(point, closest, cellId, subId, dist2)
 
-        if component_type == "proximal":
-            axial_offset = A_extra * min(0.0, np.sin(waves * base_theta))
-            distance = fabric_start_distance + axial_offset
-        elif component_type == "distal" and length != 0:
-            axial_offset = A_extra * max(0.0, np.sin(waves * base_theta))
-            distance = fabric_start_distance + length + axial_offset
-        else:
-            continue
+    cell = polydata.GetCell(cellId)
+    points = cell.GetPoints()
 
-        point, (tangent, normal, binormal) = sample_centerline_by_distance(centerline, arc_lengths, frames, distance)
-        position = point + (r + wire_radius) * (np.cos(theta_rot) * normal + np.sin(theta_rot) * binormal)
-        pts.InsertNextPoint(position)
+    point_0 = np.array(points.GetPoint(0))
+    point_1 = np.array(points.GetPoint(1))
+    point_2 = np.array(points.GetPoint(2))
 
-        if i > 0:
-            l = vtk.vtkLine()
-            l.GetPointIds().SetId(0, i-1)
-            l.GetPointIds().SetId(1, i)
-            lines.InsertNextCell(l)
+    normal = np.cross(point_1 - point_0, point_2 - point_0)
+    normal /= np.linalg.norm(normal)
 
-    if pts.GetNumberOfPoints() > 0:
-        poly = vtk.vtkPolyData()
-        poly.SetPoints(pts)
-        poly.SetLines(lines)
-        tube = vtk.vtkTubeFilter()
-        tube.SetInputData(poly)
-        tube.SetRadius(wire_radius)
-        tube.SetNumberOfSides(18)
-        tube.CappingOff()
-        tube.Update()
-        append.AddInputData(tube.GetOutput())
-
-    # ---------------------------------------------------------
-    # 2. Body Struts
-    # ---------------------------------------------------------
-    for distance_peak in distance_peaks:
-        points = vtk.vtkPoints()
-        lines_body = vtk.vtkCellArray()
+    return np.array(closest), normal
 
-        for i in range(samples + 1):
-            base_theta = 2 * np.pi * i / samples
-            theta_rot = base_theta + phase_rad
-            axial_offset = A * (1 + np.sin(waves * base_theta)) / 2.0
-            distance = distance_peak + axial_offset
+def smooth_centerline(centerline, smoothing=0.5, samples=500):
+    centerline = np.asarray(centerline)
+    spline, _ = splprep(centerline.T, s=smoothing, k=3)
+    sample_positions = np.linspace(0, 1, samples)
+    smoothed_centerline = splev(sample_positions, spline)
 
-            point, (tangent, normal, binormal) = sample_centerline_by_distance(centerline, arc_lengths, frames, distance)
-            position = point + (r + wire_radius) * (np.cos(theta_rot) * normal + np.sin(theta_rot) * binormal)
-            points.InsertNextPoint(position)
+    return np.vstack(smoothed_centerline).T
 
-            if i > 0:
-                l = vtk.vtkLine()
-                l.GetPointIds().SetId(0, i-1)
-                l.GetPointIds().SetId(1, i)
-                lines_body.InsertNextCell(l)
-
-        poly = vtk.vtkPolyData()
-        poly.SetPoints(points)
-        poly.SetLines(lines_body)
-        
-        tube = vtk.vtkTubeFilter()
-        tube.SetInputData(poly)
-        tube.SetRadius(wire_radius)
-        tube.SetNumberOfSides(18)
-        tube.CappingOff()
-        tube.Update()
-        append.AddInputData(tube.GetOutput())
+def sample_centerline_by_distance(centerline, arc_lengths, frames, query_distance):
 
-    append.Update()
-    return append.GetOutput()
+    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)
+
+    point = (1 - tangent) * point_0 + tangent * point_1
+
+    tangent_0, normal_0, binormal_0 = frames[idx-1]
+    tangent_1, normal_1, binormal_1 = 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
+
+    tangent_vector /= np.linalg.norm(tangent_vector)
+    normal_vector /= np.linalg.norm(normal_vector)
+    binormal_vector /= np.linalg.norm(binormal_vector)
+
+    return point, (tangent_vector, normal_vector, binormal_vector)
 
 # Fabric Component
 def create_fabric(centerline, frames, start_position, length, diameter, start_amplitude_pp=2.0, waves=8, phase_angle=0.0):
@@ -275,6 +155,131 @@ def create_fabric(centerline, frames, start_position, length, diameter, start_am
     
 # Struts Component
 
+def create_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, phase_angle=0.0):
+
+    arc_lengths = arc_length(centerline)
+    phase_rad = np.radians(phase_angle)
+    append = vtk.vtkAppendPolyData()
+
+    A = peak_to_peak_amplitude              # 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:]
+    elif(component_type =="extension" and (2 < len(distance_peaks) <= 4)):
+      distance_peaks = distance_peaks[2:]
+
+# fixation strut 
+    A_extra = peak_to_peak_amplitude / 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
+
+    elif component_type == "distal":
+      point, (tangent, normal, binormal) = sample_centerline_by_distance(centerline, arc_lengths, frames, fabric_start_distance + length - extra_top)
+
+    else:
+        point, (tangent, normal, binormal) = sample_centerline_by_distance(centerline, arc_lengths, frames, extra_top)
+
+    pts = vtk.vtkPoints()
+    lines = vtk.vtkCellArray()
+
+    for i in range(samples + 1):
+        base_theta = 2*np.pi*i/samples
+        theta_rot = base_theta + phase_rad
+
+        # keep only negative or positive half of sine
+        if component_type == "proximal":
+              axial_offset = A_extra * min(0.0, np.sin(waves * base_theta))
+              raw_position = (point
+                  + (r + wire_radius) * (np.cos(theta_rot)*normal + np.sin(theta_rot)*binormal)
+                  + axial_offset * tangent)
+        elif (component_type == "distal" and length != 0):
+              axial_offset = A_extra * max(0.0, np.sin(waves * base_theta))
+
+              raw_position = (point
+                        + (r + wire_radius) * (np.cos(theta_rot)*normal + np.sin(theta_rot)*binormal)
+                        + axial_offset * tangent)
+        else:
+          continue
+
+        pts.InsertNextPoint(raw_position)
+
+        if i > 0:
+            l = vtk.vtkLine()
+            l.GetPointIds().SetId(0, i-1)
+            l.GetPointIds().SetId(1, i)
+            lines.InsertNextCell(l)
+
+    poly = vtk.vtkPolyData()
+    poly.SetPoints(pts)
+    poly.SetLines(lines)
+
+    tube = vtk.vtkTubeFilter()
+    tube.SetInputData(poly)
+    tube.SetRadius(wire_radius)
+    tube.SetNumberOfSides(18)
+    tube.CappingOff()
+    tube.Update()
+
+    append.AddInputData(tube.GetOutput())
+    
+    # normal struts along the fabric body
+
+    cell_locator = vtk.vtkCellLocator()
+    cell_locator.SetDataSet(main_fabric)
+    cell_locator.BuildLocator()
+
+    points = vtk.vtkPoints()
+    lines = vtk.vtkCellArray()
+
+    for distance_peak in distance_peaks:
+
+        local_ids = []  # store point ids for this ring
+
+        for i in range(samples + 1):
+            base_theta = 2*np.pi*i/samples
+            theta_rot = base_theta + phase_rad
+            axial_offset = A * (1 + np.sin(waves * base_theta)) / 2.0
+
+            distance = distance_peak + axial_offset
+            point, (tangent, normal, binormal) = sample_centerline_by_distance(centerline, arc_lengths, frames, distance)
+
+            raw_position = (point + (r + wire_radius) * (np.cos(theta_rot)*normal + np.sin(theta_rot)*binormal))
+            surface_point, surface_normal = project_point_to_surface(raw_position, cell_locator, main_fabric)
+            position = surface_point + wire_radius * surface_normal
+           
+            point_id = points.InsertNextPoint(position)
+            local_ids.append(point_id)
+
+            if i > 0:
+                l = vtk.vtkLine()
+                l.GetPointIds().SetId(0, local_ids[i-1])
+                l.GetPointIds().SetId(1, local_ids[i])
+                lines.InsertNextCell(l)
+
+    poly = vtk.vtkPolyData()
+    poly.SetPoints(points)
+    poly.SetLines(lines)
+    tube = vtk.vtkTubeFilter()
+    tube.SetInputData(poly)
+    tube.SetRadius(wire_radius)
+    tube.SetNumberOfSides(18)
+    tube.CappingOff()
+    tube.Update()
+    append.AddInputData(tube.GetOutput())
+    append.Update()
+
+    return append.GetOutput()
 
 def create_stent_component(centerline, frames, fabric_start_distance, length, diameter, fabric, peak_to_peak=15.0, peak_spacing=7.0, component_type="proximal", phase_angle=0.0):