|
| 1 | +import numpy as np |
| 2 | +import trimesh |
| 3 | + |
| 4 | + |
| 5 | +def mesh_from_centerline( |
| 6 | + verts: np.ndarray, |
| 7 | + radii: np.ndarray, |
| 8 | + radial_segs=16, |
| 9 | + cap_segs=8, |
| 10 | + endcaps=True, |
| 11 | + is_loop=False, |
| 12 | + smooth_joints=False, |
| 13 | + joint_segs=4 |
| 14 | +) -> trimesh.Trimesh: |
| 15 | + """ |
| 16 | + Build a tube mesh with rounded ends around a polyline (rod centerline). |
| 17 | + This implementation is inspired by robust methods used in tools like Polyscope, |
| 18 | + creating clean mitered or rounded joins without resampling the centerline. |
| 19 | +
|
| 20 | + Parameters |
| 21 | + ---------- |
| 22 | + verts : (N,3) ndarray |
| 23 | + Sequence of 3D points along the rod centerline. |
| 24 | + radii : (N,) ndarray |
| 25 | + Radii at each vertex. |
| 26 | + radial_segs : int |
| 27 | + Number of segments around the tube's circumference. |
| 28 | + cap_segs : int |
| 29 | + Number of segments for the hemispherical end caps (from base to pole). |
| 30 | + endcaps : bool |
| 31 | + If True, close the ends with hemispherical caps. Ignored if is_loop is True. |
| 32 | + is_loop : bool |
| 33 | + If True, connect the ends to form a closed loop (toroid-like). |
| 34 | + smooth_joints : bool |
| 35 | + If True, creates smooth, spherical joins at vertices instead of sharp mitered joins. |
| 36 | + joint_segs : int |
| 37 | + Number of segments to use for each smoothed joint sphere. |
| 38 | +
|
| 39 | + Returns |
| 40 | + ------- |
| 41 | + mesh : trimesh.Trimesh |
| 42 | + """ |
| 43 | + verts = np.asarray(verts, dtype=float) |
| 44 | + radii = np.asarray(radii, dtype=float) |
| 45 | + N = len(verts) |
| 46 | + |
| 47 | + if N < 2: raise ValueError("Need at least 2 vertices for a rod") |
| 48 | + if is_loop and N < 3: raise ValueError("Need at least 3 vertices for a loop") |
| 49 | + if verts.shape[0] != radii.shape[0]: raise ValueError("verts and radii must have the same length") |
| 50 | + |
| 51 | + tangents = [] |
| 52 | + if is_loop: |
| 53 | + for i in range(N): |
| 54 | + tangents.append(verts[(i + 1) % N] - verts[i]) |
| 55 | + else: |
| 56 | + for i in range(N - 1): |
| 57 | + tangents.append(verts[i + 1] - verts[i]) |
| 58 | + tangents = np.array([t / np.linalg.norm(t) if np.linalg.norm(t) > 1e-9 else t for t in tangents]) |
| 59 | + |
| 60 | + # Use parallel transport to create a twist-minimizing frame |
| 61 | + basis_list = [] |
| 62 | + prev_normal = None |
| 63 | + for i in range(len(tangents)): |
| 64 | + tangent_norm = tangents[i] |
| 65 | + if i == 0: |
| 66 | + helper = np.array([0, 0, 1]) if abs(tangent_norm[2]) < 0.9 else np.array([0, 1, 0]) |
| 67 | + normal = np.cross(tangent_norm, helper) |
| 68 | + if np.linalg.norm(normal) > 1e-9: normal /= np.linalg.norm(normal) |
| 69 | + else: |
| 70 | + prev_tangent = tangents[i - 1] |
| 71 | + v = prev_tangent + tangent_norm |
| 72 | + v_dot_v = np.dot(v, v) |
| 73 | + if v_dot_v > 1e-8: |
| 74 | + reflection_vec = 2 * np.dot(prev_normal, v) / v_dot_v |
| 75 | + normal = prev_normal - v * reflection_vec |
| 76 | + else: |
| 77 | + axis = np.cross(prev_tangent, prev_normal) |
| 78 | + normal = np.cos(np.pi) * prev_normal + np.sin(np.pi) * np.cross(axis, prev_normal) |
| 79 | + |
| 80 | + binormal = np.cross(tangent_norm, normal) |
| 81 | + basis_list.append((normal, binormal)) |
| 82 | + prev_normal = normal |
| 83 | + |
| 84 | + V_list, F_list = [], [] |
| 85 | + UV_list = [] |
| 86 | + rings_in = [] |
| 87 | + rings_out = [] |
| 88 | + |
| 89 | + for i in range(len(tangents)): |
| 90 | + p_start = verts[i] |
| 91 | + p_end = verts[(i + 1) % N] |
| 92 | + r_start, r_end = radii[i], radii[(i + 1) % N] |
| 93 | + normal, binormal = basis_list[i] |
| 94 | + |
| 95 | + ring_start_indices, ring_end_indices = [], [] |
| 96 | + base_idx = len(V_list) |
| 97 | + for j in range(radial_segs): |
| 98 | + theta = 2 * np.pi * j / radial_segs |
| 99 | + offset = np.cos(theta) * normal + np.sin(theta) * binormal |
| 100 | + V_list.append(p_start + r_start * offset) |
| 101 | + V_list.append(p_end + r_end * offset) |
| 102 | + ring_start_indices.append(base_idx + 2*j) |
| 103 | + ring_end_indices.append(base_idx + 2*j + 1) |
| 104 | + |
| 105 | + u = j / (radial_segs // 2) |
| 106 | + u = 2 - u if u > 1 else u |
| 107 | + # v_start = i / max(1, len(tangents) - 1) if not is_loop else i / len(tangents) |
| 108 | + # v_end = (i + 1) / max(1, len(tangents) - 1) if not is_loop else (i + 1) / len(tangents) |
| 109 | + v_start = i |
| 110 | + v_end = i + 1 |
| 111 | + UV_list.append([u, v_start]) |
| 112 | + UV_list.append([u, v_end]) |
| 113 | + |
| 114 | + rings_in.append(ring_start_indices) |
| 115 | + rings_out.append(ring_end_indices) |
| 116 | + |
| 117 | + for j in range(radial_segs): |
| 118 | + a = ring_start_indices[j] |
| 119 | + b = ring_start_indices[(j + 1) % radial_segs] |
| 120 | + c = ring_end_indices[j] |
| 121 | + d = ring_end_indices[(j + 1) % radial_segs] |
| 122 | + F_list.extend([[a, b, c], [d, c, b]]) |
| 123 | + |
| 124 | + V_array = np.array(V_list) |
| 125 | + UV_array = np.array(UV_list) |
| 126 | + new_verts_for_joins = [] |
| 127 | + new_faces_for_joins = [] |
| 128 | + new_uvs_for_joins = [] |
| 129 | + |
| 130 | + joint_indices = range(N) if is_loop else range(1, N - 1) |
| 131 | + for i in joint_indices: |
| 132 | + idx_in = (i - 1 + N) % N |
| 133 | + idx_out = i |
| 134 | + |
| 135 | + ring_v_indices_in = rings_out[idx_in] |
| 136 | + ring_v_indices_out = rings_in[idx_out] |
| 137 | + |
| 138 | + if smooth_joints and joint_segs > 0: |
| 139 | + center, radius = verts[i], radii[i] |
| 140 | + |
| 141 | + p_in_vectors = V_array[ring_v_indices_in] - center |
| 142 | + p_in_vectors /= np.linalg.norm(p_in_vectors, axis=1, keepdims=True) |
| 143 | + p_out_vectors = V_array[ring_v_indices_out] - center |
| 144 | + p_out_vectors /= np.linalg.norm(p_out_vectors, axis=1, keepdims=True) |
| 145 | + |
| 146 | + V_array[ring_v_indices_in] = center + radius * p_in_vectors |
| 147 | + V_array[ring_v_indices_out] = center + radius * p_out_vectors |
| 148 | + |
| 149 | + omega = np.arccos(np.clip(np.dot(p_in_vectors[0], p_out_vectors[0]), -1, 1)) |
| 150 | + |
| 151 | + if omega < 1e-6 or np.isnan(omega) or abs(omega-np.pi) < 1e-6: |
| 152 | + for j in range(radial_segs): # Fallback: connect rings directly |
| 153 | + a,b,c,d = ring_v_indices_in[j], ring_v_indices_in[(j+1)%radial_segs], ring_v_indices_out[j], ring_v_indices_out[(j+1)%radial_segs] |
| 154 | + new_faces_for_joins.extend([[a, b, c], [d, c, b]]) |
| 155 | + continue |
| 156 | + |
| 157 | + sin_omega = np.sin(omega) |
| 158 | + all_rings_indices = [ring_v_indices_in] |
| 159 | + base_v_idx = len(V_array) + len(new_verts_for_joins) |
| 160 | + |
| 161 | + for k in range(1, joint_segs): |
| 162 | + t = float(k) / joint_segs |
| 163 | + current_ring_indices = [] |
| 164 | + for j in range(radial_segs): |
| 165 | + slerp_vec = (np.sin((1-t)*omega)/sin_omega) * p_in_vectors[j] + (np.sin(t*omega)/sin_omega) * p_out_vectors[j] |
| 166 | + new_verts_for_joins.append(center + radius * slerp_vec) |
| 167 | + current_ring_indices.append(base_v_idx) |
| 168 | + base_v_idx += 1 |
| 169 | + |
| 170 | + u = j / (radial_segs // 2) |
| 171 | + u = 2 - u if u > 1 else u |
| 172 | + # v = i / max(1, len(tangents) - 1) if not is_loop else i / len(tangents) |
| 173 | + v = i |
| 174 | + new_uvs_for_joins.append([u, v]) |
| 175 | + |
| 176 | + all_rings_indices.append(current_ring_indices) |
| 177 | + all_rings_indices.append(ring_v_indices_out) |
| 178 | + |
| 179 | + for k in range(len(all_rings_indices) - 1): |
| 180 | + ring_A, ring_B = all_rings_indices[k], all_rings_indices[k+1] |
| 181 | + for j in range(radial_segs): |
| 182 | + a,b,c,d = ring_A[j], ring_A[(j+1)%radial_segs], ring_B[j], ring_B[(j+1)%radial_segs] |
| 183 | + new_faces_for_joins.extend([[a, b, c], [d, c, b]]) |
| 184 | + else: # Miter join |
| 185 | + t_in, t_out = -tangents[idx_in], tangents[idx_out] |
| 186 | + miter_normal = t_in + t_out |
| 187 | + if np.linalg.norm(miter_normal) < 1e-8: continue |
| 188 | + miter_normal /= np.linalg.norm(miter_normal) |
| 189 | + for v_idx in ring_v_indices_in: |
| 190 | + p = V_array[v_idx] |
| 191 | + dist = np.dot(p - verts[i], miter_normal) |
| 192 | + V_array[v_idx] = p - dist * miter_normal |
| 193 | + for j in range(radial_segs): V_array[ring_v_indices_out[j]] = V_array[ring_v_indices_in[j]] |
| 194 | + |
| 195 | + if new_verts_for_joins: |
| 196 | + V_array = np.vstack([V_array, np.array(new_verts_for_joins)]) |
| 197 | + V_list = V_array.tolist() |
| 198 | + F_list.extend(new_faces_for_joins) |
| 199 | + if new_uvs_for_joins: |
| 200 | + UV_array = np.vstack([UV_array, np.array(new_uvs_for_joins)]) |
| 201 | + UV_list = UV_array.tolist() |
| 202 | + |
| 203 | + if not is_loop and endcaps and cap_segs > 0: |
| 204 | + for cap_type in ["start", "end"]: |
| 205 | + is_start = cap_type == "start" |
| 206 | + center, radius = (verts[0], radii[0]) if is_start else (verts[-1], radii[-1]) |
| 207 | + tangent = -tangents[0] if is_start else tangents[-1] |
| 208 | + normal, binormal = basis_list[0] if is_start else basis_list[-1] |
| 209 | + prev_ring_indices = rings_in[0] if is_start else rings_out[-1] |
| 210 | + |
| 211 | + for k in range(1, cap_segs + 1): |
| 212 | + alpha, is_pole = k * (np.pi/2) / cap_segs, k == cap_segs |
| 213 | + ring_radius, displacement = radius * np.cos(alpha), radius * np.sin(alpha) |
| 214 | + ring_center = center + displacement * tangent |
| 215 | + current_ring_indices = [] |
| 216 | + if not is_pole: |
| 217 | + for j in range(radial_segs): |
| 218 | + theta = 2*np.pi * j / radial_segs |
| 219 | + offset = np.cos(theta)*normal + np.sin(theta)*binormal |
| 220 | + V_list.append(ring_center + ring_radius*offset) |
| 221 | + current_ring_indices.append(len(V_list)-1) |
| 222 | + |
| 223 | + u = j / radial_segs |
| 224 | + v = ring_radius / radius |
| 225 | + UV_list.append([u, v]) |
| 226 | + else: |
| 227 | + V_list.append(ring_center) |
| 228 | + current_ring_indices = [len(V_list)-1] * radial_segs |
| 229 | + UV_list.append([0.0, 0.0]) # Center of endcap |
| 230 | + |
| 231 | + for j in range(radial_segs): |
| 232 | + a,b = prev_ring_indices[j], prev_ring_indices[(j+1)%radial_segs] |
| 233 | + c,d = current_ring_indices[j], current_ring_indices[(j+1)%radial_segs] |
| 234 | + if not is_pole: |
| 235 | + faces = [[a,c,b], [d,b,c]] if is_start else [[a,b,c], [d,c,b]] |
| 236 | + F_list.extend(faces) |
| 237 | + else: |
| 238 | + F_list.append([b,c,a] if is_start else [a,b,c]) |
| 239 | + prev_ring_indices = current_ring_indices |
| 240 | + |
| 241 | + mesh = trimesh.Trimesh( |
| 242 | + vertices=np.array(V_list), |
| 243 | + faces=np.array(F_list, dtype=int), |
| 244 | + visual=trimesh.visual.TextureVisuals(uv=np.array(UV_list)), |
| 245 | + process=True, |
| 246 | + ) |
| 247 | + return mesh |
0 commit comments