辅助资料：正文课程已提供中英双语。本页保存当前源码或示例检查记录；代码标识符与原始结果不作改写。先读概念与推导，再用本页逐项对照。
Supporting material: the lessons now include Chinese and English explanations. This page retains source excerpts or verification records; identifiers and original numerical results remain unchanged. Learn the concepts and derivations first, then compare them with these records.
把公式对照到实际函数
以下内容从当前 Python 文件原样摘录，并保留行号。阅读注释时仍要核对循环的实际顺序：逐点更新未处理的相邻段可能暂时失真；“硬约束”应在完整扫描和最终结果中检查。
这些是课程生成时的源码快照；后续修改实现后需要重新生成。ShapeIK 当前最终恢复使用最短弧，广义 Helix 当前使用一维二分与数值积分。
src/dhik/solvers/fabrik.py · FABRIK._iterate · 第 272 行
272 @staticmethod
273 def _iterate(
274 positions: np.ndarray,
275 base: np.ndarray,
276 target: np.ndarray,
277 bones: np.ndarray,
278 rest_directions: np.ndarray,
279 max_iterations: int,
280 tolerance: float,
281 ) -> tuple:
282 """执行 FABRIK 前后向迭代，返回 ``(实际迭代次数, 是否收敛)``。
283
284 ``positions`` 会被**原地更新**。每一步都严格把每段投影回 ``bones`` 给定的长度，
285 因此任意中间状态都不改变骨长；段方向未定义时用该段的静止骨向兜底。
286 """
287 num_bones = int(bones.shape[0])
288 reached = False
289 iterations = 0
290 for iteration in range(1, max_iterations + 1):
291 iterations = iteration
292 # 前向：末端拉到目标，再从末端往根部投影。
293 positions[num_bones] = target
294 for step in range(num_bones - 1, -1, -1):
295 segment = positions[step] - positions[step + 1]
296 unit = _unit_or_fallback(segment, rest_directions[step])
297 positions[step] = positions[step + 1] + unit * bones[step]
298 # 后向：根部复位，再从根部往末端投影。
299 positions[0] = base
300 for step in range(1, num_bones + 1):
301 segment = positions[step] - positions[step - 1]
302 unit = _unit_or_fallback(segment, rest_directions[step - 1])
303 positions[step] = positions[step - 1] + unit * bones[step - 1]
304 if float(np.linalg.norm(positions[num_bones] - target)) <= tolerance:
305 reached = True
306 break
307 return iterations, reached
src/dhik/solvers/shape.py · ShapeIK.solve · 第 228 行
228 def solve(
229 self,
230 skeleton: Skeleton,
231 target: np.ndarray,
232 pole: Optional[np.ndarray] = None,
233 end_effector: Optional[JointRef] = None,
234 root_transform: Optional[Transform] = None,
235 parameters: Optional[SolverParameters] = None,
236 ) -> SolveResult:
237 """见 :meth:`IKSolver.solve`。
238
239 与其它求解器的差异：除了末端位置，本求解器还决定整条链的**形状**
240 （由 ``family`` 选形状族、``pole`` 选弯曲平面与弯曲侧、``α`` 由端点和
241 骨长自动解出）。不可达 / 退化 / 未收敛都通过 ``status`` 报告。
242 """
243 start = self.stopwatch_start()
244 values = self.resolve_parameters(parameters)
245 index, point, pole_vector = self.validate_common_inputs(skeleton, target, pole, end_effector)
246 try:
247 spec = _SolverSpec.from_parameters(values)
248 except _ParameterError as error:
249 return self.build_failure_result(
250 skeleton,
251 point,
252 SolveStatus.INVALID_INPUT,
253 message=str(error),
254 root_transform=root_transform,
255 end_effector=index,
256 )
257
258 chain = skeleton.ancestors(index)
259 num_bones = len(chain) - 1
260 if num_bones < 1:
261 raise ValueError(
262 "ShapeIK 需要至少 1 段骨的骨链，但末端就是骨架根关节（没有可驱动的骨）；"
263 f"请指定 root 之外的末端关节，当前末端为 {skeleton.joint_names[index]!r}。"
264 )
265 bones = np.array([skeleton.bone_length(chain[j]) for j in range(1, len(chain))], dtype=np.float64)
266 total_length = float(bones.sum())
267 radius_min, radius_max = reachable_bounds(skeleton, index)
268
269 root = Transform.identity() if root_transform is None else root_transform
270 base = root.apply(skeleton.rest_position_of(chain[0]))
271 offset = point - base
272 distance = float(np.linalg.norm(offset))
273 scale = max(1.0, total_length)
274 # 弦方向：端点重合时方向未定义，用确定性垂直方向兜底（不产生 NaN）。
275 chord_hat = (
276 offset / distance if distance > _DIRECTION_EPS else orthogonal_unit_vector(np.zeros(3))
277 )
278 bend, bend_note = _bend_direction(chord_hat, pole_vector, 1e-9 * scale)
279 # 形状族在拿到弯曲方向之后才构造：plane_hint 就是"曲线朝哪侧凸出"，
280 # 于是换 pole 就换弯曲平面/弯曲侧（形状参数 α 与平面无关，见 shapes.py 约定）。
281 family = spec.build_family(bend)
282 classification = classify_target(skeleton, point, end_effector=index, root_transform=root)
283
284 messages: List[str] = []
285 if bend_note:
286 messages.append(bend_note)
287 if classification is Reachability.NEAR_STRAIGHT:
288 messages.append("目标位于外边界邻域（近完全伸直），按可达处理但形状参数已逼近奇异")
289 elif classification is Reachability.NEAR_FOLDED:
290 messages.append("目标位于内边界邻域（近完全折叠），按可达处理但形状参数已逼近奇异")
291
292 # --- 形状参数求解 + 等弧长采样初始化 --------------------------------- #
293 parameter: Optional[Any] = None
294 initial: Optional[np.ndarray] = None
295 degenerate: Optional[ShapeDegenerateError] = None
296 try:
297 parameter = family.solve_parameter(
298 chord=distance, total_length=total_length, lengths=bones, **spec.control
299 )
300 cumulative = np.concatenate(([0.0], np.cumsum(bones)))
301 initial, polyline = _shape_sampling(family, base, point, parameter, cumulative)
302 except ShapeDegenerateError as error:
303 degenerate = error
304 else:
305 bulge = _max_bulge(initial, base, chord_hat)
306 limit = _MIN_SHAPE_BULGE_RATIO * scale
307 if bulge <= limit:
308 degenerate = ShapeDegenerateError(
309 f"形状参数 {family.parameter_to_metadata(parameter)!r} 解出的曲线与弦的最大垂距 "
310 f"{bulge:.3e} ≤ {limit:.3e}（浮点精度下不可分辨，曲线已退化成直线）",
311 "减小形状控制量（正弦族可减小波数 k）或改用其它形状族",
312 )
313
314 # --- 离散骨长修正（硬约束）+ 关节旋转恢复 ---------------------------- #
315 correction: Optional[CorrectionResult] = None
316 fallback: Optional[str] = None
317 initial_residual: Optional[float] = None
318 if degenerate is None:
319 assert parameter is not None and initial is not None # 非退化分支必然已解出参数与采样
320 initial_residual = _max_segment_residual(initial, bones)
321 projector = _CurveProjector(polyline)
322 correction = correct_to_bone_lengths(
323 points=initial,
324 lengths=bones,
325 A=base,
326 B=point,
327 max_iter=spec.max_iterations,
328 tol=spec.tolerance,
329 shape_projector=projector,
330 attraction=spec.attraction,
331 )
332 points = correction.points
333 iterations = correction.iterations
334 status = SolveStatus.SUCCESS if correction.converged else SolveStatus.MAX_ITERATIONS
335 if not correction.converged:
336 messages.append(
337 f"修正器在 {correction.iterations} 轮内未收敛到容差 {spec.tolerance:.3g}"
338 f"（端点偏移 {correction.endpoint_error:.3e}、段长残差 "
339 f"{correction.max_segment_residual:.3e}）"
340 )
341 else:
342 messages.append(f"形状参数退化：{degenerate.reason}；建议回退：{degenerate.fallback_hint}")
343 points, iterations, status, fallback, correction, initial_residual = _fallback_configuration(
344 classification=classification,
345 base=base,
346 point=point,
347 chord_hat=chord_hat,
348 bones=bones,
349 bend=bend,
350 spec=spec,
351 )
352 if fallback == "straight":
353 messages.append("已回退为沿弦方向的伸直姿态（骨长仍严格保持）")
354 elif fallback == "folded":
355 messages.append("已回退为折叠极限姿态（最长骨朝目标方向、其余骨折回，末端落在最小可达半径上）")
356 else:
357 messages.append("已放松形状先验为直线并求可行姿态（骨长仍严格保持）")
358
359 local_rotations = self.recover_local_rotations(skeleton, chain, points, root)
360 metadata = _build_metadata(
361 spec=spec,
362 family=family,
363 classification=classification,
364 parameter=parameter,
365 num_bones=num_bones,
366 bones=bones,
367 total_length=total_length,
368 radius_min=radius_min,
369 radius_max=radius_max,
370 distance=distance,
371 bend=bend,
372 pole_vector=pole_vector,
373 points=points,
374 base=base,
375 point=point,
376 chord_hat=chord_hat,
377 correction=correction,
378 initial_residual=initial_residual,
379 degenerate=degenerate,
380 fallback=fallback,
381 )
382 elapsed = self.stopwatch_elapsed(start)
383 return self.build_result(
384 skeleton=skeleton,
385 local_rotations=local_rotations,
386 target=point,
387 status=status,
388 iterations=iterations,
389 elapsed_time=elapsed,
390 root_transform=root,
391 end_effector=index,
392 message="；".join(messages),
393 metadata=metadata,
394 )
src/dhik/shapes.py · ArcShape.solve_parameter · 第 615 行
615 def solve_parameter(self, chord: float, total_length: float, **control: Any) -> float:
616 """解圆心角 ``θ``：见类文档的单调性与二分区间说明。"""
617 distance, total = self._validate_inputs(chord, total_length, dict(control))
618 ratio = total / distance
619 upper = float(np.pi) * (1.0 - 1e-12)
620 if _arc_length_ratio(upper) < ratio:
621 raise ShapeDegenerateError(
622 f"圆心角需要逼近 2π 才能满足 S/d = {ratio:.6g}（超过可表示范围 ~1e12）",
623 "回退为多圈形状或退化为接近闭合的圆环姿态",
624 )
625 half_angle = _bisect_increasing(_arc_length_ratio, 0.0, upper, ratio)
626 return 2.0 * half_angle
src/dhik/shapes.py · CatenaryShape.solve_parameter · 第 713 行
713 def solve_parameter(self, chord: float, total_length: float, **control: Any) -> float:
714 """解 ``u = d/(2a)`` 后返回 ``a = d/(2u)``：见类文档的单调性与二分区间说明。"""
715 distance, total = self._validate_inputs(chord, total_length, dict(control))
716 ratio = total / distance
717 upper = 1.0
718 while _catenary_length_ratio(upper) < ratio:
719 upper *= 2.0
720 if upper > 700.0:
721 raise ShapeDegenerateError(
722 f"悬链线需要 u = d/(2a) > 700 才能满足 S/d = {ratio:.6g}（sinh 溢出）",
723 "回退为沿弦方向的伸直姿态，或改用圆弧/正弦族表达极端弯曲",
724 )
725 u = _bisect_increasing(_catenary_length_ratio, 0.0, upper, ratio)
726 return distance / (2.0 * u)
src/dhik/shapes.py · SinusoidShape._points_at_fractions · 第 874 行
874 def _points_at_fractions(
875 self,
876 A: np.ndarray,
877 B: np.ndarray,
878 parameter: float,
879 fractions: np.ndarray,
880 ) -> np.ndarray:
881 """按弧长比例取点：对 ``s(v) = Q(v)/(kπ)`` 二分反解出 ``v``，再代回 ``p(t)``。
882
883 ``Q(v) = ∫₀^v √(d² + c²cos²u) du`` 用"整周期解析 + 余项复合 Gauss-Legendre"
884 计算（见 :func:`_sinusoid_arc_integral`），避免 ``c ≫ d`` 时 ``|cos|`` 拐点
885 带来的求积误差。
886 """
887 chord_hat, _, distance = self._chord(A, B)
888 amplitude = _as_scalar(parameter, "parameter（正弦幅值 α）")
889 plane = self._plane_direction(A, B)
890 fractions = np.asarray(fractions, dtype=np.float64).reshape(-1)
891 k = self.k
892 c = amplitude * k * np.pi
893 v_max = float(k) * np.pi
894 total = self._arc_length_of(distance, c)
895 if not total > 0.0:
896 raise ShapeDegenerateError(
897 f"正弦曲线弧长非正（S = {total:.6g}），无法按弧长采样",
898 "检查幅值与端点是否合法，或改用其它形状族",
899 )
900
901 def cumulative(v: float) -> float:
902 return _sinusoid_arc_integral(distance, c, v) / v_max
903
904 v = np.empty_like(fractions)
905 for index, fraction in enumerate(fractions):
906 if fraction <= 0.0:
907 v[index] = 0.0
908 elif fraction >= 1.0:
909 v[index] = v_max
910 else:
911 v[index] = _bisect_increasing(cumulative, 0.0, v_max, float(fraction) * total, iterations=64)
912 t = v / v_max
913 offset = (B - A)[None, :] * t[:, None]
914 bulge = (amplitude * np.sin(np.pi * k * t))[:, None] * plane[None, :]
915 return A[None, :] + offset + bulge
src/dhik/shapes.py · correct_to_bone_lengths · 第 1031 行
1031 def correct_to_bone_lengths(
1032 points: np.ndarray,
1033 lengths: np.ndarray,
1034 A: np.ndarray,
1035 B: np.ndarray,
1036 max_iter: int = 200,
1037 tol: float = 1e-12,
1038 shape_projector: Optional[Callable[[np.ndarray], np.ndarray]] = None,
1039 attraction: float = 0.5,
1040 ) -> CorrectionResult:
1041 """把等弧长采样得到的关节点修正到**严格满足离散骨长约束**。
1042
1043 为什么需要它（课程点名的致命陷阱）
1044 ----------------------------------
1045 "曲线弧长不等于相邻关节点间的直线距离"：等弧长采样只保证相邻点之间的**弧长**
1046 等于骨长，而约束要求的是**直线距离** ``‖p_{i+1} − p_i‖ = l_i``。曲线上两点间的
1047 弦永远短于它们之间的弧，所以等弧长采样天然偏短（弯得越厉害偏得越多）。
1048 本函数以采样结果为初始化，迭代拉回硬约束。
1049
1050 算法（Gauss-Seidel 交替投影，FABRIK 式）
1051 ----------------------------------------
1052 每一轮做以下三步，返回前**最后一步一定是纯硬约束扫描**（吸引步只在
1053 非最后一轮执行，且每轮都以一次硬投影收尾，因此返回构型永远严格满足骨长）：
1054
1055 1. **前向**（保 ``A`` 固定）：``p_0 ← A``，逐段把 ``p_{i+1}`` 投到以 ``p_i`` 为心、
1056 半径 ``l_i`` 的球面上——即沿当前方向重新量出 ``l_i``；
1057 2. **形状吸引**（可选，仅当给出 ``shape_projector`` 且 ``attraction > 0``）：
1058 每个内点朝曲线最近点混一步 ``p ← (1 − w)·p + w·proj(p)``（形状是**软约束**，
1059 这一步会破坏骨长，由下一步立刻修回）；该步只在**形状偏置阶段**执行
1060 （前 ``min(20, max_iter // 2)`` 轮），让软约束决定解落在哪一支，
1061 随后交给纯硬约束扫描收敛；
1062 3. **后向**（保 ``B`` 固定）：``p_N ← B``，从末端反扫，同样按 ``l_i`` 重摆每段。
1063
1064 由于每步都把段长写成半径投影，**任意中间构型的段长都精确等于 ``l_i``**
1065 （只受浮点舍入影响），这正是"骨长是硬约束、绝不静默改变骨长"的体现；
1066 收敛判据是端点误差与段长残差都进入 ``tol``。
1067
1068 形状偏置阶段为什么要短（≤ 20 轮）而且必须收尾
1069 ----------------------------------------------
1070 若每轮都吸引，软约束会把构型持续推向曲线；当曲线本身接近直线时，构型会贴到
1071 **共线流形**上——而共线方向正是投影扫描的不动流形（投影只沿线重排，不产生
1072 横向弯曲），此时硬约束的收敛率退化到接近 1：实测端点误差长期停在 ``1e-2``
1073 量级（须几万轮才收敛）。因此偏置只做一个短的前置阶段，之后全是纯硬约束；
1074 这样既保住了形状偏置带来的解分支（形状偏差实测下降 3 成以上），又能把
1075 端点误差压回机器精度。
1076
1077 参数:
1078 points: 形状 ``(N, 3)`` 的初始关节点（等弧长采样结果），``N = len(lengths) + 1``。
1079 lengths: 形状 ``(N−1,)`` 的骨长，必须全为正。
1080 A: 形状 ``(3,)`` 的起点（硬约束，逐轮钉住）。
1081 B: 形状 ``(3,)`` 的终点（硬约束）。
1082 max_iter: 最大扫描轮数（每轮 = 吸引 + 前向 + 后向）。
1083 tol: 收敛容差（端点误差与段长残差同时满足）。
1084 shape_projector: 可选的可调用对象，接收 ``(M, 3)`` 的内点数组并返回其
1085 在曲线上的最近点 ``(M, 3)``；为 ``None`` 时不吸引、偏差字段为 ``None``。
1086 attraction: 吸引步的混合系数 ``w ∈ [0, 1]``（0 等价于不吸引），
1087 只在前 ``min(20, max_iter // 2)`` 轮生效（形状偏置阶段，见下方说明）。
1088
1089 返回:
1090 :class:`CorrectionResult`（含修正前后形状偏差对比）。
1091
1092 异常:
1093 ValueError: 输入形状/取值非法（骨长非正、端点形状错、``attraction`` 越界、
1094 ``shape_projector`` 返回形状不符等）。**不可达（``‖B − A‖`` 小于骨链
1095 最小可达半径）不抛异常**：如实通过 ``converged=False`` 与
1096 ``endpoint_error > 0`` 报告，骨长仍然精确。
1097 """
1098 start = _as_point(A, "A")
1099 end = _as_point(B, "B")
1100 current = np.asarray(points, dtype=np.float64)
1101 if current.ndim != 2 or current.shape[1] != 3:
1102 raise ValueError(f"points 形状必须为 (N, 3)，收到 {current.shape}")
1103 if not bool(np.all(np.isfinite(current))):
1104 raise ValueError("points 含 NaN 或无穷值")
1105 bones = np.asarray(lengths, dtype=np.float64).reshape(-1)
1106 if bones.size == 0:
1107 raise ValueError("lengths 不能为空（至少 1 段骨）")
1108 if current.shape[0] != bones.size + 1:
1109 raise ValueError(
1110 f"points 的关节点数应等于骨数 + 1：收到 {current.shape[0]} 个点、{bones.size} 根骨"
1111 )
1112 if not bool(np.all(np.isfinite(bones))) or bool(np.any(bones <= 0.0)):
1113 raise ValueError(f"lengths 必须是有限正数，收到 {bones.tolist()}")
1114 if not np.isfinite(float(bones.sum())):
1115 # 防止骨长量级过大时逐段累加溢出成 Inf（宁可明确报错，也不静默产出 Inf）。
1116 raise ValueError(f"lengths 之和溢出（骨长量级过大），收到 {bones.tolist()}")
1117 rounds = int(max_iter)
1118 if rounds < 1:
1119 raise ValueError(f"max_iter 必须 ≥ 1，收到 {max_iter!r}")
1120 tolerance = _as_scalar(tol, "tol")
1121 if tolerance <= 0.0:
1122 raise ValueError(f"tol 必须为正，收到 {tol!r}")
1123 weight = _as_scalar(attraction, "attraction")
1124 if not 0.0 <= weight <= 1.0:
1125 raise ValueError(f"attraction 必须落在 [0, 1]，收到 {attraction!r}")
1126 if shape_projector is not None and not callable(shape_projector):
1127 raise TypeError("shape_projector 必须可调用（接收 (M, 3) 返回 (M, 3)）或为 None")
1128
1129 num_bones = int(bones.size)
1130 current = current.copy()
1131 chord = end - start
1132 chord_hat = _unit_or(chord, _deterministic_perpendicular(chord))
1133 deviation_before = (
1134 _shape_deviation(current, shape_projector) if shape_projector is not None else None
1135 )
1136
1137 def forward_sweep() -> None:
1138 """前向：``A`` 固定，逐段把下一点投到以当前点为心、半径 ``l_i`` 的球面上。"""
1139 current[0] = start
1140 direction = chord_hat
1141 for index in range(num_bones):
1142 direction = _unit_or(current[index + 1] - current[index], direction)
1143 current[index + 1] = current[index] + direction * bones[index]
1144
1145 def backward_sweep() -> None:
1146 """后向：``B`` 固定，从末端反扫恢复端点（FABRIK 式）。"""
1147 current[-1] = end
1148 direction = -chord_hat
1149 for index in range(num_bones - 1, -1, -1):
1150 direction = _unit_or(current[index] - current[index + 1], direction)
1151 current[index] = current[index + 1] + direction * bones[index]
1152
1153 def attract() -> None:
1154 """形状吸引步（软约束）：内点朝曲线最近点混合 ``weight``。"""
1155 if num_bones < 2:
1156 return
1157 interior = current[1:-1]
1158 projected = _project_points(shape_projector, interior)
1159 current[1:-1] = interior + weight * (projected - interior)
1160
1161 # 形状偏置阶段取前 min(20, max_iter//2) 轮：软约束负责把构型拉进"贴合形状"的
1162 # 那一支，其后（含最后一轮）全是纯硬约束扫描，把端点误差收敛到 0。
1163 # 偏置阶段不能太长：实测把构型拉到近共线流形后，硬约束扫描的收敛率会退化到
1164 # 接近 1（共线方向是投影的不动流形），端点误差会长期停在 1e-2 量级上不去。
1165 bias_rounds = min(20, max(1, rounds // 2))
1166 iterations = 0
1167 converged = False
1168 for round_index in range(1, rounds + 1):
1169 iterations = round_index
1170 # 最后一轮不吸引：返回前必有一次纯硬约束扫描收尾。
1171 if (
1172 shape_projector is not None
1173 and weight > 0.0
1174 and round_index <= bias_rounds
1175 and round_index < rounds
1176 ):
1177 attract()
1178 forward_sweep()
1179 backward_sweep()
1180 endpoint_error = _endpoint_error(current, start, end)
1181 residual = _max_segment_residual(current, bones)
1182 if endpoint_error <= tolerance and residual <= tolerance:
1183 converged = True
1184 break
1185
1186 endpoint_error = _endpoint_error(current, start, end)
1187 residual = _max_segment_residual(current, bones)
1188 deviation_after = (
1189 _shape_deviation(current, shape_projector) if shape_projector is not None else None
1190 )
1191 return CorrectionResult(
1192 points=current,
1193 iterations=iterations,
1194 converged=converged,
1195 endpoint_error=endpoint_error,
1196 max_segment_residual=residual,
1197 shape_deviation_before=deviation_before,
1198 shape_deviation_after=deviation_after,
1199 )
src/dhik/solvers/shape.py · _CurveProjector.__call__ · 第 556 行
556 def __call__(self, points: np.ndarray) -> np.ndarray:
557 query = np.asarray(points, dtype=np.float64)
558 if query.size == 0:
559 return query.copy()
560 starts = self._polyline[:-1]
561 segments = self._polyline[1:] - starts
562 length_squared = np.einsum("ij,ij->i", segments, segments)
563 relative = query[:, None, :] - starts[None, :, :]
564 projection = np.einsum("mij,ij->mi", relative, segments)
565 safe = np.where(length_squared > 0.0, length_squared, 1.0)
566 ratio = np.clip(projection / safe[None, :], 0.0, 1.0)
567 candidates = starts[None, :, :] + ratio[:, :, None] * segments[None, :, :]
568 distances = np.linalg.norm(query[:, None, :] - candidates, axis=2)
569 nearest = np.argmin(distances, axis=1)
570 return candidates[np.arange(query.shape[0]), nearest]
src/dhik/shapes_helix.py · HelixShape.solve_parameter · 第 251 行
251 def solve_parameter(self, chord: float, total_length: float, **control: Any) -> HelixParameter:
252 """由弦长 ``d`` 与目标弧长 ``L`` 解出标准或广义螺旋参数。
253
254 标准版：螺距显式解，半径解析解。广义版：显式解平均螺距，
255 再在非负半径域内二分求平均半径；还须满足 |Δr| <= d、L >= S_min。
256
257 参数:
258 chord: 弦长 ``d = ‖B − A‖``。
259 total_length: 目标弧长 ``L``（取骨链总长）。
260 **control: 认 ``turns``（覆盖实例圈数，须与实例一致）、``sign``、
261 ``lengths``（骨长序列，用于内不可达判定）；``taper_radius`` /
262 ``taper_pitch`` 覆盖本次半径和螺距的线性变化量。
263
264 返回:
265 :class:`HelixParameter`。
266
267 异常:
268 ShapeDegenerateError: ``d ≈ 0`` / ``L ≤ d`` / ``L ≈ d`` / ``d < r_min``。
269 ValueError: 圈数非正整数、非有限 taper、控制量与实例不一致等。
270 """
271 distance, total = self._validate_inputs(chord, total_length, dict(control))
272 turns = self._resolve_turns(control)
273 sign = self._resolve_sign(control)
274 dr = _finite(control.get("taper_radius", self.taper_radius), "taper_radius")
275 dp = _finite(control.get("taper_pitch", self.taper_pitch), "taper_pitch")
276 if dr != 0.0 or dp != 0.0:
277 return self._solve_general(distance, total, turns, sign, dr, dp)
278 pbar = sign * distance / float(turns)
279 inner = (total / float(turns)) ** 2 - pbar * pbar
280 if inner <= 0.0:
281 # L > d 已由 _validate_inputs 保证，这里只是数值兜底。
282 raise ShapeDegenerateError(
283 f"标准 helix 在 N = {turns} 下需要 L/N > |p̄| = {abs(pbar):.6g}，"
284 f"当前 L/N = {total / float(turns):.6g} 不满足",
285 "增大骨长或减小弦长；或降低圈数 N",
286 )
287 radius = float(np.sqrt(inner) / (2.0 * np.pi))
288 if not np.isfinite(radius):
289 raise ShapeDegenerateError(
290 f"解出的半径非有限（L = {total:.6g}、d = {distance:.6g}、N = {turns}）",
291 "检查输入量级，或换用其它形状族",
292 )
293 return HelixParameter(r0=radius, r1=radius, p0=pbar, p1=pbar, turns=float(turns))
src/dhik/shapes_helix.py · HelixShape._solve_general · 第 448 行
448 def _solve_general(self, distance, total, turns, sign, dr, dp):
449 """先解弦长，再在 r_bar >= |dr|/2 的单调区间二分。"""
450 if abs(dr) > distance:
451 raise ShapeDegenerateError("半径变化量超过端点距离，当前螺旋形状无解", "减小 taper_radius")
452 # 避免 d**2 - dr**2 的相减损失。
453 pbar = sign * np.sqrt((distance - abs(dr)) * (distance + abs(dr))) / turns
454 p0, p1 = pbar - 0.5 * dp, pbar + 0.5 * dp
455 low = 0.5 * abs(dr)
456 def length(radius):
457 return _arc_length_closed(radius - 0.5 * dr, radius + 0.5 * dr, p0, p1, turns)
458 minimum = length(low)
459 tolerance = 1e-12 * max(1.0, total)
460 if minimum > total + tolerance:
461 raise ShapeDegenerateError("该锥度的最短曲线仍超过骨链总长", "减小 taper 或增大骨链总长")
462 if abs(minimum - total) <= tolerance:
463 radius = low
464 else:
465 high = max(total, low + total)
466 if length(high) < total:
467 raise ShapeDegenerateError("无法建立有限求根区间", "检查参数量级")
468 for _ in range(60):
469 mid = 0.5 * (low + high)
470 if length(mid) < total:
471 low = mid
472 else:
473 high = mid
474 if high - low <= 1e-13 * max(1.0, total):
475 break
476 radius = 0.5 * (low + high)
477 return HelixParameter(max(0.0, radius - 0.5 * dr), max(0.0, radius + 0.5 * dr),
478 p0, p1, float(turns))
src/dhik/shapes_helix.py · _arc_at · 第 157 行
157 def _arc_at(r0: float, r1: float, p0: float, p1: float, turns: float, u: float) -> float:
158 """累计弧长：标准版解析，广义版复合 Gauss-Legendre 积分。"""
159 if r0 == r1 and p0 == p1:
160 return float(np.hypot(2.0 * np.pi * turns * r0, turns * p0) * u)
161 t = float(u) * _ARC_U
162 dr, dp = r1 - r0, p1 - p0
163 speed = np.sqrt(dr * dr + (2.0 * np.pi * turns * (r0 + dr * t)) ** 2
164 + (turns * (p0 + dp * t)) ** 2)
165 return float(u) * float(np.dot(_ARC_W, speed))
src/dhik/solvers/base.py · IKSolver.recover_local_rotations · 第 386 行
386 @staticmethod
387 def recover_local_rotations(
388 skeleton: Skeleton,
389 chain_indices: Sequence[int],
390 chain_positions: np.ndarray,
391 root_transform: Optional[Transform] = None,
392 ) -> np.ndarray:
393 """由骨链关节的**世界位置**恢复**局部关节旋转向量**（位置 → 旋转）。
394
395 项目要求"求解结果必须恢复为局部关节旋转，再通过 FK 重建骨架"，本函数就是
396 这条要求的公共实现：迭代型/形状型求解器算出关节点位置后调用它，解析型
397 求解器推导出关节点位置后同样调用它，从而两类方法口径一致。
398
399 推导（FK 递推 ``p_child = p_parent + R_parent · t_child^rest``）：
400
401 对链上相邻的父子关节 ``(i, j)``，有 ``R_i · t_j^rest = p_j − p_i``。
402 已知 ``R_i = R_parent(i) · R(θ_i)``，于是
403
404 .. math::
405
406 R(θ_i) \\, t_j^{rest} = \\big(R_{parent(i)}\\big)^{T} (p_j - p_i)
407
408 右边只依赖已求出的父朝向与两点位置，问题化为"把一个已知方向旋到另一个
409 已知方向"，用 :func:`minimal_rotation_vector`（最短弧）解出 ``θ_i``。
410 沿链从根到末端一次扫描即可，**骨长自动保持**（旋转不改变 ``t_j^rest`` 的模长）。
411
412 扭转自由度口径：每个关节只约束"骨轴指向"，绕骨轴的自转是自由的；这里统一
413 取**最短弧**（不引入额外轴向扭转），使同一位置序列永远给出同一组关节角。
414 不在链上的关节（含末端关节自身）局部旋转保持为零。
415
416 参数:
417 skeleton: 骨架。
418 chain_indices: 骨链关节下标序列（根 → 末端，拓扑序，必须从骨架根关节开始）。
419 chain_positions: 形状 ``(len(chain_indices), 3)`` 的各关节世界位置。
420 root_transform: 施加在根关节上的世界变换（决定链根的父朝向）。
421
422 返回:
423 形状 ``(N, 3)`` 的局部关节旋转向量。
424
425 异常:
426 ValueError: 下标序列为空、不以骨架根关节开头、长度或位置形状不匹配，
427 或相邻关节点重合（方向未定义）。
428 """
429 indices = [int(i) for i in chain_indices]
430 if not indices:
431 raise ValueError("chain_indices 不能为空")
432 if indices[0] != skeleton.root_index:
433 raise ValueError(
434 "chain_indices 必须从骨架根关节开始（根 → 末端），"
435 f"收到起点下标 {indices[0]}，根关节下标为 {skeleton.root_index}"
436 )
437 positions = np.asarray(chain_positions, dtype=np.float64)
438 if positions.shape != (len(indices), 3):
439 raise ValueError(
440 f"chain_positions 形状必须为 ({len(indices)}, 3)，收到 {positions.shape}"
441 )
442 if not np.all(np.isfinite(positions)):
443 raise ValueError("chain_positions 包含 NaN 或无穷值")
444
445 root = Transform.identity() if root_transform is None else root_transform
446 local_rotations = np.zeros((skeleton.num_joints, 3), dtype=np.float64)
447 # accumulated = 链根关节**父级**的世界朝向（之后逐关节右乘自身旋转）。
448 accumulated = root.rotation.copy()
449 for offset in range(len(indices) - 1):
450 parent = indices[offset]
451 child = indices[offset + 1]
452 bone_vector = positions[offset + 1] - positions[offset]
453 rest_offset = skeleton.rest_offset_of(child)
454 # accumulated^T · bone_vector 就是 bone_vector 在父关节坐标系下的表达。
455 target_in_parent = accumulated.T @ bone_vector
456 # 先用最短弧给出"把该段骨转过去"的相对旋转，再经 rotation.py 的
457 # "旋转矩阵 → 旋转向量"（Shepperd 分支，θ ∈ [0, π] 规范形式）取回关节角，
458 # 保证扭转自由度口径与 FK/旋转模块完全一致。
459 relative_matrix = rotation_vector_to_matrix(
460 minimal_rotation_vector(rest_offset, target_in_parent)
461 )
462 local_rotations[parent] = matrix_to_rotation_vector(relative_matrix)
463 accumulated = accumulated @ relative_matrix
464 return local_rotations