Developable Surfaces

1import numpy as np
2
3def unroll_developable_strip(strip_triangles):
4    """
5    Isometrically unfolds a 3D developable triangle strip onto the 2D XY plane
6    by rotating each successive triangle around its shared ruling edge via Rodrigues' formula.
7    """
8    def rodrigues_rotate(vec, axis, angle):
9        axis = axis / np.linalg.norm(axis)
10        return (vec * np.cos(angle) +
11                np.cross(axis, vec) * np.sin(angle) +
12                axis * np.dot(axis, vec) * (1.0 - np.cos(angle)))
13
14    unrolled = [np.asarray(strip_triangles[0], dtype=float)]
15    for k in range(1, len(strip_triangles)):
16        prev_tri = unrolled[-1]
17        curr_tri = np.asarray(strip_triangles[k], dtype=float)
18        # Shared ruling edge between (prev_tri[1], prev_tri[2]) and (curr_tri[0], curr_tri[1])
19        n_prev = np.cross(prev_tri[1] - prev_tri[0], prev_tri[2] - prev_tri[0])
20        n_curr = np.cross(curr_tri[1] - curr_tri[0], curr_tri[2] - curr_tri[0])
21        n_prev /= np.linalg.norm(n_prev)
22        n_curr /= np.linalg.norm(n_curr)
23        hinge_axis = prev_tri[2] - prev_tri[1]
24        dihedral = np.arctan2(np.dot(np.cross(n_curr, n_prev), hinge_axis / np.linalg.norm(hinge_axis)),
25                              np.dot(n_curr, n_prev))
26        tip_vec = rodrigues_rotate(curr_tri[2] - curr_tri[0], hinge_axis, dihedral)
27        unrolled.append(np.array([prev_tri[1], prev_tri[2], prev_tri[1] + tip_vec]))
28    return unrolled

1import numpy as np
2
3def compute_geodesic_curvature(curve_points, surface_normals):
4    """
5    Computes discrete geodesic curvature kappa_g along a surface curve.
6    Decomposes curve acceleration into normal curvature kappa_n (along N)
7    and geodesic curvature kappa_g (along tangent binormal B_g = N x T).
8    A timber lamella is geodesic iff kappa_g == 0 everywhere.
9    """
10    pts = np.asarray(curve_points, dtype=float)
11    normals = np.asarray(surface_normals, dtype=float)
12    kappa_g = np.zeros(len(pts) - 2)
13
14    for i in range(1, len(pts) - 1):
15        e_prev = pts[i] - pts[i - 1]
16        e_next = pts[i + 1] - pts[i]
17        ds = 0.5 * (np.linalg.norm(e_prev) + np.linalg.norm(e_next))
18        t_prev = e_prev / np.linalg.norm(e_prev)
19        t_next = e_next / np.linalg.norm(e_next)
20
21        # Discrete curvature vector k * N_curve = dT / ds
22        curvature_vec = (t_next - t_prev) / max(ds, 1e-8)
23        tangent = (t_prev + t_next) / np.linalg.norm(t_prev + t_next)
24        surf_normal = normals[i] / np.linalg.norm(normals[i])
25        tangent_normal = np.cross(surf_normal, tangent)
26        kappa_g[i - 1] = np.dot(curvature_vec, tangent_normal)
27    return kappa_g