Skip to main content

animsmith_core/
metrics.rs

1//! Locomotion clip metrics: loop-seam ratio, gait phase, root-motion
2//! speed, and sampled root trajectory. The loop-seam, gait, and speed
3//! metrics were ported from a production game pipeline's reference
4//! implementation
5//! (verified there against Blender pose-matrix FK to <0.01×) — the
6//! algorithms are kept semantically identical so the numbers reproduce.
7
8use crate::model::{Clip, Document, Property, Track};
9use crate::profile::{ResolvedRoles, Role};
10use crate::sample::{PoseGrid, sample_clip};
11use glam::{DQuat, DVec3, Quat, Vec3};
12use serde::{Deserialize, Serialize};
13use std::cell::RefCell;
14use std::collections::BTreeMap;
15use std::rc::Rc;
16
17/// Below this per-frame foot move (metres), a clip has no real stride
18/// (idle / block / stationary action) and the seam ratio would be a
19/// divide-by-noise, so no ratio is reported.
20pub const MIN_STRIDE_STEP_M: f64 = 0.02;
21
22/// Half-turn tolerance used to refuse direction-ambiguous adjacent yaw steps
23/// and to preserve the unwrapped sign when canonicalizing multi-step endpoint
24/// results near ±180 degrees.
25pub const ROOT_YAW_HALF_TURN_AMBIGUITY_DEG: f64 = 1.0e-4;
26
27/// Lazily sampled metric pose grids for one document.
28///
29/// The check, measurement, and report pipelines all judge the same
30/// uniform metric grid. Sharing this owner lets callers run checks and
31/// then emit measurements or reports without sampling the same clip
32/// twice.
33///
34/// The cache uses `Rc` and `RefCell`, so it is intentionally neither
35/// `Send` nor `Sync`. Create one owner per document on each worker thread,
36/// then share it by reference among consumers on that thread.
37#[derive(Debug)]
38pub struct MetricGrids<'a> {
39    doc: &'a Document,
40    grids: RefCell<BTreeMap<usize, Rc<PoseGrid>>>,
41}
42
43impl<'a> MetricGrids<'a> {
44    /// Create a lazy metric-grid cache for `doc`.
45    pub fn new(doc: &'a Document) -> Self {
46        Self {
47            doc,
48            grids: RefCell::new(BTreeMap::new()),
49        }
50    }
51
52    /// The document these grids sample.
53    pub fn document(&self) -> &'a Document {
54        self.doc
55    }
56
57    /// The metric pose grid for clip `clip_index`, computed once and
58    /// shared. Returns `None` for an out-of-range index, non-positive
59    /// duration, or fewer than three keys on the longest track.
60    pub fn grid(&self, clip_index: usize) -> Option<Rc<PoseGrid>> {
61        let clip = self.doc.clips.get(clip_index)?;
62        let frames = metric_frame_count(clip)?;
63        Some(
64            self.grids
65                .borrow_mut()
66                .entry(clip_index)
67                .or_insert_with(|| Rc::new(sample_clip(&self.doc.skeleton, clip, frames)))
68                .clone(),
69        )
70    }
71}
72
73/// Foot-cycle metrics for one sampled clip.
74#[derive(Debug, Clone, PartialEq)]
75#[non_exhaustive]
76pub struct FootCycleMetrics {
77    /// Wrap discontinuity of the feet (relative to hips) over the max of
78    /// the two seam-adjacent in-clip steps. ≈1.0 for a clean cyclic
79    /// loop; well above 1 for a seam pop. `None` when the clip has no
80    /// real stride (see [`Self::has_real_stride`]), or when a real
81    /// stride exists but `seam / neighbour_step` is not finite. The
82    /// whole-clip position finiteness gate above only bounds individual
83    /// coordinates, not the `f32` squared-length arithmetic used to turn
84    /// two positions into a distance; a per-axis delta near `f32::MAX`
85    /// overflows that squaring to infinity even though every input
86    /// position was finite. That is the only known route to this case —
87    /// it requires magnitudes far outside any real animation.
88    pub loop_seam_ratio: Option<f64>,
89    /// Whether the seam-adjacent neighbour step met the configured
90    /// minimum stride threshold, i.e. whether the clip has a real
91    /// stride to normalize the seam against. `false` means
92    /// [`Self::loop_seam_ratio`]'s absence is an expected "no subject"
93    /// (a planted/idle clip), not a derivation failure.
94    pub has_real_stride: bool,
95    /// Cycle position `[0,1)` of the trough of the fundamental harmonic of the
96    /// left-minus-right foot-height signal — a stride-phase anchor encoding
97    /// handedness + cycle alignment. `None` when a side is missing, the
98    /// sampled signal has exact zero peak-to-peak swing, or the harmonic fit
99    /// fails.
100    pub gait_phase: Option<f64>,
101    /// Peak-to-peak swing of the L−R foot-height signal (metres); near
102    /// zero means no detectable alternation and the phase is noise.
103    pub lr_amplitude_m: f64,
104}
105
106/// Why a clip does or does not carry a stride anchor.
107///
108/// This is the vocabulary the gait checks report coverage in, and the one a
109/// presentation reads to say why a member of a declared group has no anchor.
110#[derive(Debug, Clone, Copy, PartialEq)]
111#[non_exhaustive]
112pub enum GaitPhaseOutcome {
113    /// One side's foot roles did not resolve, so there is no left-minus-right
114    /// signal to fit a phase to.
115    MissingBilateralFootRoles,
116    /// The left-minus-right signal has exactly zero peak-to-peak swing, so
117    /// its phase would be noise rather than a stride.
118    NoFootHeightSwing,
119    /// The stride anchor, in cycle fraction `[0, 1)`.
120    Measured(f64),
121    /// The harmonic fit yielded no finite phase.
122    Unavailable,
123}
124
125impl FootCycleMetrics {
126    /// Why this clip does or does not carry a stride anchor.
127    pub fn gait_phase_outcome(&self, roles: &ResolvedRoles) -> GaitPhaseOutcome {
128        GaitPhaseOutcome::classify(
129            self.gait_phase,
130            self.lr_amplitude_m,
131            bilateral_foot_roles(roles),
132        )
133    }
134}
135
136/// Whether both sides resolve a foot or toe role, which is what a
137/// left-minus-right signal needs.
138fn bilateral_foot_roles(roles: &ResolvedRoles) -> bool {
139    let side = |foot, toe| roles.get(foot).is_some() || roles.get(toe).is_some();
140    side(Role::LeftFoot, Role::LeftToe) && side(Role::RightFoot, Role::RightToe)
141}
142
143/// Whether the rig resolves the roles any gait measurement needs: hips plus
144/// at least one foot.
145///
146/// This is the prerequisite the gait checks report as a coverage gap when it
147/// fails; [`gait_member_phase`] applies it, so a consumer reads the outcome
148/// rather than the predicate.
149pub(crate) fn gait_roles_resolved(roles: &ResolvedRoles) -> bool {
150    let has_foot = [
151        Role::LeftFoot,
152        Role::LeftToe,
153        Role::RightFoot,
154        Role::RightToe,
155    ]
156    .iter()
157    .any(|&role| roles.get(role).is_some());
158    roles.get(Role::Hips).is_some() && has_foot
159}
160
161impl GaitPhaseOutcome {
162    fn classify(gait_phase: Option<f64>, lr_amplitude_m: f64, has_bilateral_roles: bool) -> Self {
163        match gait_phase {
164            _ if !has_bilateral_roles => Self::MissingBilateralFootRoles,
165            _ if lr_amplitude_m == 0.0 => Self::NoFootHeightSwing,
166            Some(phase) => Self::Measured(phase),
167            None => Self::Unavailable,
168        }
169    }
170}
171
172/// Model-space loop-continuity measurements for one skeleton bone.
173#[derive(Debug, Clone, PartialEq)]
174#[non_exhaustive]
175pub struct BoneLoopContinuityMetrics {
176    /// Last-sample to first-sample model-space position distance (metres).
177    pub position_delta_m: f64,
178    /// Shortest-path model-space rotation difference (degrees).
179    pub rotation_delta_deg: f64,
180    /// Difference between the model-space linear velocities immediately
181    /// before and after the wrap (metres per second).
182    pub seam_velocity_delta_mps: f64,
183    /// Difference between the model-space angular velocities immediately
184    /// before and after the wrap (degrees per second).
185    pub seam_angular_velocity_delta_degps: f64,
186}
187
188/// Sampled model-space translation and yaw for a selected Root/Hips bone.
189#[derive(Debug, Clone, PartialEq)]
190#[non_exhaustive]
191pub struct RootTrajectoryMetrics {
192    /// Translation facts when every selected-bone position is finite.
193    pub translation: Option<RootTranslationMetrics>,
194    /// Yaw facts when a fixed, deterministic horizontal heading witness remains
195    /// usable across the complete sampled trajectory.
196    pub yaw: Option<RootYawMetrics>,
197}
198
199/// Sampled model-space translation facts for a selected Root/Hips bone.
200#[derive(Debug, Clone, Copy, PartialEq)]
201#[non_exhaustive]
202pub struct RootTranslationMetrics {
203    /// Endpoint displacement along canonical model-space +X, in metres.
204    pub horizontal_displacement_x_m: f64,
205    /// Endpoint displacement along canonical model-space +Z, in metres.
206    pub horizontal_displacement_z_m: f64,
207    /// Sum of sampled model-space XZ step lengths, in metres.
208    pub horizontal_travel_m: f64,
209    /// Signed endpoint displacement along canonical model-space +Y, in metres.
210    pub vertical_displacement_m: f64,
211    /// Minimum signed +Y displacement from the initial sample, in metres.
212    pub vertical_min_displacement_m: f64,
213    /// Maximum signed +Y displacement from the initial sample, in metres.
214    pub vertical_max_displacement_m: f64,
215}
216
217/// Signed sampled yaw facts for a selected Root/Hips bone.
218#[derive(Debug, Clone, Copy, PartialEq)]
219#[non_exhaustive]
220pub struct RootYawMetrics {
221    /// Fixed local basis axis used as the horizontal heading witness.
222    pub heading_axis: RootYawHeadingAxis,
223    /// Shortest signed endpoint-equivalent yaw in `[-180, 180]` degrees.
224    /// An exact half turn retains the sign of [`Self::unwrapped_yaw_deg`].
225    pub net_yaw_deg: f64,
226    /// Signed first-to-last heading change after deterministic wrap crossing
227    /// unwrapping. Unlike endpoint orientation alone, a sampled full turn is
228    /// retained as approximately `+360` or `-360` degrees.
229    pub unwrapped_yaw_deg: f64,
230    /// Sum of absolute sampled unwrapped heading steps. This retains reversing
231    /// yaw motion that cancels in [`Self::unwrapped_yaw_deg`].
232    pub yaw_travel_deg: f64,
233}
234
235/// One local orientation witness retained for a complete model-space yaw
236/// measurement. Order is policy: conventional `+Z` wins an exact tie,
237/// followed by the common Z-up-source `+Y` convention and finally `+X`.
238#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
239#[serde(rename_all = "snake_case")]
240#[non_exhaustive]
241pub enum RootYawHeadingAxis {
242    /// Positive local Z axis.
243    PositiveZ,
244    /// Positive local Y axis.
245    PositiveY,
246    /// Positive local X axis.
247    PositiveX,
248}
249
250impl RootYawHeadingAxis {
251    /// Stable display label for the selected local heading witness.
252    pub const fn label(self) -> &'static str {
253        match self {
254            Self::PositiveZ => "+Z",
255            Self::PositiveY => "+Y",
256            Self::PositiveX => "+X",
257        }
258    }
259}
260
261/// Model-space horizontal `(x, z)` projection of one local unit basis axis.
262pub(crate) fn horizontal_heading(rotation: DQuat, axis: RootYawHeadingAxis) -> (f64, f64) {
263    let local_axis = match axis {
264        RootYawHeadingAxis::PositiveZ => DVec3::Z,
265        RootYawHeadingAxis::PositiveY => DVec3::Y,
266        RootYawHeadingAxis::PositiveX => DVec3::X,
267    };
268    let heading = rotation.mul_vec3(local_axis);
269    (heading.x, heading.z)
270}
271
272/// Select the best-conditioned local heading witness at sample zero. Exact
273/// ties retain the declared `+Z`, `+Y`, `+X` priority because replacement is
274/// strict rather than `>=`.
275pub(crate) fn select_horizontal_heading_axis(rotation: DQuat) -> RootYawHeadingAxis {
276    let mut selected = RootYawHeadingAxis::PositiveZ;
277    let (x, z) = horizontal_heading(rotation, selected);
278    let mut best_length = x.hypot(z);
279    for candidate in [RootYawHeadingAxis::PositiveY, RootYawHeadingAxis::PositiveX] {
280        let (x, z) = horizontal_heading(rotation, candidate);
281        let length = x.hypot(z);
282        if length > best_length {
283            selected = candidate;
284            best_length = length;
285        }
286    }
287    selected
288}
289
290/// Convert signed unwrapped yaw to its canonical endpoint-equivalent value.
291/// Endpoint-equivalent results within
292/// [`ROOT_YAW_HALF_TURN_AMBIGUITY_DEG`] of a half turn retain the sign of the
293/// unwrapped result so analytic `+180` and `-180` trajectories remain
294/// distinguishable despite binary32 quaternion roundoff.
295pub(crate) fn canonical_net_yaw_deg(unwrapped_yaw_deg: f64) -> f64 {
296    let mut net = (unwrapped_yaw_deg + 180.0).rem_euclid(360.0) - 180.0;
297    if (net.abs() - 180.0).abs() <= ROOT_YAW_HALF_TURN_AMBIGUITY_DEG {
298        net = 180.0f64.copysign(unwrapped_yaw_deg);
299    }
300    net
301}
302
303/// Measure model-space translation and sampled yaw for `bone`.
304///
305/// AnimSmith's metric domain is right-handed, +Y-up metres; horizontal
306/// displacement is therefore the signed X/Z endpoint vector and vertical
307/// evidence is signed +Y displacement from sample zero. The existing uniform
308/// [`PoseGrid`] is authoritative. Yaw chooses at sample zero whichever local
309/// `+Z`, `+Y`, or `+X` axis has the greatest horizontal projection, retains
310/// that witness for every sample, and unwraps crossings at ±180 degrees.
311/// Positive yaw increases `atan2(x, z)`; for a +Z-aligned witness this rotates
312/// +Z toward +X, the positive right-handed direction around normalized +Y. A multi-step
313/// result within [`ROOT_YAW_HALF_TURN_AMBIGUITY_DEG`] of a half turn is
314/// canonicalized to signed ±180 using the unwrapped sign.
315/// Exact 180-degree adjacent steps are ambiguous and make only yaw unavailable.
316///
317/// Returns `None` when the selected bone is outside the grid or fewer than two
318/// samples exist. Translation and yaw are derived independently: non-finite
319/// positions set [`RootTrajectoryMetrics::translation`] to `None`, while
320/// rotation or heading failures set [`RootTrajectoryMetrics::yaw`] to `None`.
321pub fn root_trajectory_metrics(grid: &PoseGrid, bone: usize) -> Option<RootTrajectoryMetrics> {
322    let frames = grid.frame_count();
323    if frames < 2 || bone >= grid.bone_count() {
324        return None;
325    }
326
327    let first = grid.model_position(0, bone);
328    let mut translation_valid = first.is_finite();
329    let first_x = f64::from(first.x);
330    let first_y = f64::from(first.y);
331    let first_z = f64::from(first.z);
332    let mut last_x = first_x;
333    let mut last_y = first_y;
334    let mut last_z = first_z;
335    let mut previous_x = first_x;
336    let mut previous_z = first_z;
337    let mut horizontal_travel_m = 0.0f64;
338    let mut vertical_min_displacement_m = 0.0f64;
339    let mut vertical_max_displacement_m = 0.0f64;
340
341    let mut yaw_valid = true;
342    let mut heading_axis = None;
343    let mut first_heading_deg: Option<f64> = None;
344    let mut previous_heading_deg: Option<f64> = None;
345    let mut winding_turns = 0i64;
346    let mut yaw_travel_deg = 0.0f64;
347
348    for frame in 0..frames {
349        let position = grid.model_position(frame, bone);
350        if !position.is_finite() {
351            translation_valid = false;
352        } else if translation_valid {
353            last_x = f64::from(position.x);
354            last_y = f64::from(position.y);
355            last_z = f64::from(position.z);
356            if frame > 0 {
357                horizontal_travel_m += (last_x - previous_x).hypot(last_z - previous_z);
358            }
359            previous_x = last_x;
360            previous_z = last_z;
361            let vertical = last_y - first_y;
362            vertical_min_displacement_m = vertical_min_displacement_m.min(vertical);
363            vertical_max_displacement_m = vertical_max_displacement_m.max(vertical);
364        }
365
366        if !yaw_valid {
367            continue;
368        }
369        let rotation = grid.model_rotation(frame, bone);
370        let length_squared = rotation.length_squared();
371        if !rotation.is_finite() || !length_squared.is_finite() || length_squared == 0.0 {
372            yaw_valid = false;
373            continue;
374        }
375        let rotation = rotation.as_dquat().normalize();
376        let axis = *heading_axis.get_or_insert_with(|| select_horizontal_heading_axis(rotation));
377        let (heading_x, heading_z) = horizontal_heading(rotation, axis);
378        let horizontal_length = heading_x.hypot(heading_z);
379        if !horizontal_length.is_finite() || horizontal_length <= f64::from(f32::EPSILON) {
380            yaw_valid = false;
381            continue;
382        }
383        let heading_deg = heading_x.atan2(heading_z).to_degrees();
384        if let Some(previous) = previous_heading_deg {
385            let raw_delta = heading_deg - previous;
386            if (raw_delta.abs() - 180.0).abs() <= ROOT_YAW_HALF_TURN_AMBIGUITY_DEG {
387                yaw_valid = false;
388                continue;
389            }
390            if raw_delta > 180.0 {
391                winding_turns -= 1;
392                yaw_travel_deg += (raw_delta - 360.0).abs();
393            } else if raw_delta < -180.0 {
394                winding_turns += 1;
395                yaw_travel_deg += (raw_delta + 360.0).abs();
396            } else {
397                yaw_travel_deg += raw_delta.abs();
398            }
399        } else {
400            first_heading_deg = Some(heading_deg);
401        }
402        previous_heading_deg = Some(heading_deg);
403    }
404
405    // `PoseGrid` positions are binary32 and the grid cannot exceed
406    // `usize::MAX` frames, so finite samples cannot overflow these widened
407    // binary64 endpoint, extrema, or accumulated-travel calculations. Keep
408    // the final filter as a fail-closed boundary if either representation
409    // changes; model-space FK overflow is rejected in the loop above.
410    let translation = translation_valid
411        .then_some(RootTranslationMetrics {
412            horizontal_displacement_x_m: last_x - first_x,
413            horizontal_displacement_z_m: last_z - first_z,
414            horizontal_travel_m,
415            vertical_displacement_m: last_y - first_y,
416            vertical_min_displacement_m,
417            vertical_max_displacement_m,
418        })
419        .filter(|translation| {
420            [
421                translation.horizontal_displacement_x_m,
422                translation.horizontal_displacement_z_m,
423                translation.horizontal_travel_m,
424                translation.vertical_displacement_m,
425                translation.vertical_min_displacement_m,
426                translation.vertical_max_displacement_m,
427            ]
428            .into_iter()
429            .all(f64::is_finite)
430        });
431
432    let yaw = yaw_valid.then(|| {
433        let unwrapped_yaw_deg = previous_heading_deg.expect("non-empty pose grid")
434            - first_heading_deg.expect("non-empty pose grid")
435            + winding_turns as f64 * 360.0;
436        RootYawMetrics {
437            heading_axis: heading_axis.expect("non-empty pose grid"),
438            net_yaw_deg: canonical_net_yaw_deg(unwrapped_yaw_deg),
439            unwrapped_yaw_deg,
440            yaw_travel_deg,
441        }
442    });
443    let yaw = yaw.filter(|yaw| {
444        yaw.net_yaw_deg.is_finite()
445            && yaw.unwrapped_yaw_deg.is_finite()
446            && yaw.yaw_travel_deg.is_finite()
447    });
448
449    Some(RootTrajectoryMetrics { translation, yaw })
450}
451
452/// Return the shortest-path model-space rotation vector from `from` to `to`.
453///
454/// The left-relative step (`to * from⁻¹`) expresses the angular direction in
455/// model space. Canonicalizing the quaternion hemisphere makes the result
456/// invariant to the equivalent `q`/`-q` representation. At exactly 180
457/// degrees, where `w` cannot choose a hemisphere, the first non-zero vector
458/// component breaks the tie deterministically.
459fn shortest_path_model_rotation_vector(from: Quat, to: Quat) -> Option<Vec3> {
460    let mut step = to * from.conjugate();
461    if !step.is_finite() {
462        return None;
463    }
464
465    let [x, y, z, w] = step.to_array();
466    if w < 0.0 || (w == 0.0 && (x < 0.0 || (x == 0.0 && (y < 0.0 || (y == 0.0 && z < 0.0))))) {
467        step = -step;
468    }
469
470    let vector = step.xyz();
471    let sin_half_angle = vector.length();
472    if !sin_half_angle.is_finite() {
473        return None;
474    }
475    if sin_half_angle == 0.0 {
476        return Some(Vec3::ZERO);
477    }
478
479    let angle_rad = 2.0 * sin_half_angle.atan2(step.w);
480    let rotation_vector = vector * (angle_rad / sin_half_angle);
481    rotation_vector.is_finite().then_some(rotation_vector)
482}
483
484/// Measure C0 pose closure plus C1 linear- and angular-velocity continuity
485/// independently for every bone.
486///
487/// The grid spans `[0, duration]`, including both endpoints. C1 continuity is
488/// therefore the difference between the in-clip step entering the last sample
489/// and the in-clip step leaving frame 0. Treating the last-to-first endpoint
490/// chord as a velocity would assign zero velocity to a perfectly closed loop.
491///
492/// Returns `None` when the shared grid has fewer than three frames, has no
493/// bones, or has an unusable seam-adjacent time step. A row is `None` only
494/// when that bone's seam-adjacent model-space evidence is unusable; one bad
495/// bone never suppresses finite evidence for another bone.
496pub fn loop_continuity_metrics(grid: &PoseGrid) -> Option<Vec<Option<BoneLoopContinuityMetrics>>> {
497    let frames = grid.frame_count();
498    if frames < 3 || grid.bone_count() == 0 {
499        return None;
500    }
501
502    let first_dt = f64::from(grid.times[1] - grid.times[0]);
503    let last_dt = f64::from(grid.times[frames - 1] - grid.times[frames - 2]);
504    if !first_dt.is_finite() || !last_dt.is_finite() || first_dt <= 0.0 || last_dt <= 0.0 {
505        return None;
506    }
507
508    Some(
509        (0..grid.bone_count())
510            .map(|bone| {
511                let first = grid.model_position(0, bone);
512                let next = grid.model_position(1, bone);
513                let previous = grid.model_position(frames - 2, bone);
514                let last = grid.model_position(frames - 1, bone);
515                if [first, next, previous, last]
516                    .iter()
517                    .any(|position| !position.is_finite())
518                {
519                    return None;
520                }
521
522                let rotations = [
523                    grid.model_rotation(0, bone),
524                    grid.model_rotation(1, bone),
525                    grid.model_rotation(frames - 2, bone),
526                    grid.model_rotation(frames - 1, bone),
527                ];
528                if rotations.iter().any(|rotation| {
529                    !rotation.is_finite()
530                        || !rotation.length_squared().is_finite()
531                        || rotation.length_squared() == 0.0
532                }) {
533                    return None;
534                }
535                let [
536                    first_rotation,
537                    next_rotation,
538                    previous_rotation,
539                    last_rotation,
540                ] = rotations.map(Quat::normalize);
541                let delta = first_rotation.conjugate() * last_rotation;
542                let [x, y, z, w] = delta.to_array();
543                let sin_half_angle = Vec3::new(x, y, z).length();
544                let rotation_delta_deg =
545                    f64::from(2.0 * sin_half_angle.atan2(w.abs()).to_degrees());
546                let position_delta_m = f64::from((last - first).length());
547                let outgoing_velocity = (next - first) / first_dt as f32;
548                let incoming_velocity = (last - previous) / last_dt as f32;
549                let seam_velocity_delta_mps =
550                    f64::from((outgoing_velocity - incoming_velocity).length());
551                let outgoing_angular_velocity =
552                    shortest_path_model_rotation_vector(first_rotation, next_rotation)?
553                        / first_dt as f32;
554                let incoming_angular_velocity =
555                    shortest_path_model_rotation_vector(previous_rotation, last_rotation)?
556                        / last_dt as f32;
557                let seam_angular_velocity_delta_degps = f64::from(
558                    (outgoing_angular_velocity - incoming_angular_velocity)
559                        .length()
560                        .to_degrees(),
561                );
562
563                if !position_delta_m.is_finite()
564                    || !rotation_delta_deg.is_finite()
565                    || !seam_velocity_delta_mps.is_finite()
566                    || !seam_angular_velocity_delta_degps.is_finite()
567                {
568                    return None;
569                }
570                Some(BoneLoopContinuityMetrics {
571                    position_delta_m,
572                    rotation_delta_deg,
573                    seam_velocity_delta_mps,
574                    seam_angular_velocity_delta_degps,
575                })
576            })
577            .collect(),
578    )
579}
580
581/// Measure the foot cycle of a clip from its pose grid. Requires the
582/// Hips role and at least one foot role; returns `None` otherwise (the
583/// caller decides which typed coverage gap represents the missing metric).
584///
585/// The grid must span `[0, duration]` — the wrap pair is
586/// `(last frame, frame 0)`. Grids under 3 frames carry no cycle.
587///
588/// # Panics
589///
590/// Panics if `roles` contains bone indices outside `grid`. Role
591/// resolutions produced by this crate are tied to the same skeleton that
592/// produced the grid; embedders that hand-build roles must preserve that
593/// relationship.
594pub fn foot_cycle_metrics(
595    grid: &PoseGrid,
596    roles: &ResolvedRoles,
597    min_stride_step_m: f64,
598) -> Option<FootCycleMetrics> {
599    let GaitBones { hips, left, right } = gait_bones(grid, roles)?;
600    let feet: Vec<usize> = left.iter().chain(right.iter()).copied().collect();
601    let frames = grid.frame_count();
602    // Feet relative to hips: cancels the in-place root so we measure
603    // the leg cycle, not body travel.
604    let rel = |frame: usize, bone: usize| -> Vec3 {
605        grid.model_position(frame, bone) - grid.model_position(frame, hips)
606    };
607
608    // Loop seam: the wrap chord vs its NEIGHBOURING in-clip steps (the
609    // step into the last frame and the step out of the first) — local
610    // continuity, because stride speed varies legitimately inside a
611    // cycle and the wrap may sit at an arbitrary cycle position. A real
612    // pop is discontinuous against its immediate neighbours too.
613    let max_foot_dist = |a: usize, b: usize| -> f64 {
614        feet.iter()
615            .map(|&f| (rel(a, f) - rel(b, f)).length() as f64)
616            .fold(0.0, f64::max)
617    };
618    let seam = max_foot_dist(frames - 1, 0);
619    let step_first = max_foot_dist(1, 0);
620    let step_last = max_foot_dist(frames - 1, frames - 2);
621    let neighbour_step = step_first.max(step_last);
622    let has_real_stride = neighbour_step > 0.0 && neighbour_step >= min_stride_step_m;
623    let loop_seam_ratio = if has_real_stride {
624        let ratio = seam / neighbour_step;
625        ratio.is_finite().then_some(ratio)
626    } else {
627        None
628    };
629
630    // Gait phase: fundamental-harmonic trough of the L−R foot-height
631    // signal over one cycle (the duplicate wrap frame excluded). The
632    // difference cancels common-mode pelvis bob and encodes handedness
633    // plus a stable cycle anchor. One function forms that signal and fits
634    // it, so a caller that draws the curve draws what was measured.
635    let evidence = lr_evidence(grid, hips, &left, &right);
636    Some(FootCycleMetrics {
637        loop_seam_ratio,
638        has_real_stride,
639        gait_phase: match evidence.outcome {
640            GaitPhaseOutcome::Measured(phase) => Some(phase),
641            _ => None,
642        },
643        lr_amplitude_m: evidence.lr_amplitude_m,
644    })
645}
646
647/// The bones a gait measurement reads, once the grid is known to carry a
648/// cycle it can be read from.
649struct GaitBones {
650    hips: usize,
651    left: Vec<usize>,
652    right: Vec<usize>,
653}
654
655/// Resolve the bones a gait measurement needs from `grid` and `roles`, or
656/// `None` where the clip carries no readable foot cycle at all: under three
657/// frames, no hips, no foot on either side, or a non-finite sample anywhere
658/// in the hips-relative feet.
659///
660/// This is the single gate [`foot_cycle_metrics`] and [`gait_phase_evidence`]
661/// share, so a caller cannot see a phase the other would refuse to measure.
662fn gait_bones(grid: &PoseGrid, roles: &ResolvedRoles) -> Option<GaitBones> {
663    if grid.frame_count() < 3 {
664        return None;
665    }
666    let hips = roles.get(Role::Hips)?;
667    let side =
668        |foot, toe| -> Vec<usize> { [foot, toe].iter().filter_map(|&r| roles.get(r)).collect() };
669    let left = side(Role::LeftFoot, Role::LeftToe);
670    let right = side(Role::RightFoot, Role::RightToe);
671    let feet: Vec<usize> = left.iter().chain(right.iter()).copied().collect();
672    if feet.is_empty() {
673        return None;
674    }
675    let finite = (0..grid.frame_count()).all(|frame| {
676        let hips_position = grid.model_position(frame, hips);
677        hips_position.is_finite()
678            && feet
679                .iter()
680                .all(|&foot| (grid.model_position(frame, foot) - hips_position).is_finite())
681    });
682    finite.then_some(GaitBones { hips, left, right })
683}
684
685/// The stride-anchor evidence one clip yields.
686///
687/// The signal, the samples a cycle spans, the fit's outcome and the swing an
688/// amplitude floor is judged against all come from one measurement, so a
689/// drawn curve, a drawn anchor and a reported coverage gap describe the same
690/// thing rather than three re-derivations of it.
691#[derive(Debug, Clone, PartialEq)]
692#[non_exhaustive]
693pub struct GaitPhaseEvidence {
694    /// Left-minus-right foot height relative to the hips, in metres, one
695    /// value per sampled frame of the clip's metric grid. Empty when the rig
696    /// resolves a foot on one side only, which is a left-minus-right signal
697    /// with no right-hand term.
698    pub lr_foot_height_m: Vec<f64>,
699    /// Samples one stride cycle spans, so sample `k` of
700    /// [`Self::lr_foot_height_m`] sits at cycle position `k / cycle_samples`
701    /// and a phase `p` sits at sample `p * cycle_samples`. A grid over three
702    /// frames excludes its duplicate wrap sample from the cycle; a
703    /// three-frame grid has no duplicate to exclude.
704    pub cycle_samples: usize,
705    /// Why this clip does or does not carry a stride anchor.
706    pub outcome: GaitPhaseOutcome,
707    /// Peak-to-peak swing of the signal over one cycle (metres).
708    pub lr_amplitude_m: f64,
709}
710
711/// Measure one clip's stride-anchor evidence from its metric grid.
712///
713/// `None` where the clip carries no readable foot cycle at all — the same
714/// refusal [`foot_cycle_metrics`] makes, through the same gate.
715///
716/// # Panics
717///
718/// Panics if `roles` contains bone indices outside `grid`, exactly as
719/// [`foot_cycle_metrics`] does.
720pub fn gait_phase_evidence(grid: &PoseGrid, roles: &ResolvedRoles) -> Option<GaitPhaseEvidence> {
721    let GaitBones { hips, left, right } = gait_bones(grid, roles)?;
722    Some(lr_evidence(grid, hips, &left, &right))
723}
724
725/// Samples one stride cycle spans on a `frames`-sample metric grid.
726///
727/// A grid of more than three frames ends on a repeat of its first sample —
728/// the wrap — and that duplicate is not part of the cycle, so the cycle is
729/// `frames - 1`. A three-frame grid carries no such duplicate, and all three
730/// of its samples are the cycle. Sample `k` therefore sits at cycle position
731/// `k / cycle`, which is where a phase measured on that cycle is drawn
732/// against it.
733pub fn gait_cycle_samples(frames: usize) -> usize {
734    if frames > 3 { frames - 1 } else { frames }
735}
736
737/// Form the left-minus-right foot-height signal and fit its stride anchor.
738fn lr_evidence(grid: &PoseGrid, hips: usize, left: &[usize], right: &[usize]) -> GaitPhaseEvidence {
739    let frames = grid.frame_count();
740    let cycle_samples = gait_cycle_samples(frames);
741    if left.is_empty() || right.is_empty() {
742        return GaitPhaseEvidence {
743            lr_foot_height_m: Vec::new(),
744            cycle_samples,
745            outcome: GaitPhaseOutcome::MissingBilateralFootRoles,
746            lr_amplitude_m: 0.0,
747        };
748    }
749    let mean_height = |frame: usize, bones: &[usize]| -> f64 {
750        let hips_position = grid.model_position(frame, hips);
751        bones
752            .iter()
753            .map(|&bone| (grid.model_position(frame, bone) - hips_position).y as f64)
754            .sum::<f64>()
755            / bones.len() as f64
756    };
757    let lr_foot_height_m: Vec<f64> = (0..frames)
758        .map(|frame| mean_height(frame, left) - mean_height(frame, right))
759        .collect();
760    // The amplitude and the fit read one cycle; the series keeps every
761    // sampled frame, because a caller drawing the curve draws the clip.
762    let cycle = &lr_foot_height_m[..cycle_samples];
763    let max = cycle.iter().copied().fold(f64::MIN, f64::max);
764    let min = cycle.iter().copied().fold(f64::MAX, f64::min);
765    let lr_amplitude_m = max - min;
766    let phase = (lr_amplitude_m > 0.0)
767        .then(|| fundamental_trough_phase(cycle))
768        .flatten();
769    GaitPhaseEvidence {
770        lr_foot_height_m,
771        cycle_samples,
772        outcome: GaitPhaseOutcome::classify(phase, lr_amplitude_m, true),
773        lr_amplitude_m,
774    }
775}
776
777/// How one declared member of a gait group stands, classified once for every
778/// consumer.
779///
780/// The `gait-group` check turns these into findings and coverage gaps; a
781/// report draws the measured ones and names the rest. Deriving them twice is
782/// what lets a picture claim a member was measured where the check recorded
783/// a gap.
784#[derive(Debug, Clone, Copy, PartialEq)]
785#[non_exhaustive]
786pub enum GaitMemberPhase {
787    /// The document holds no clip of this member's name.
788    Absent,
789    /// The rig did not resolve hips together with a foot role on at least
790    /// one side — either prerequisite being absent is enough — so no member
791    /// of any group has a phase to measure.
792    RolesUnresolved,
793    /// The clip is present but carries no readable foot cycle.
794    NoFootCycle,
795    /// The foot cycle was read and yielded no anchor, with the reason.
796    NoAnchor(GaitPhaseOutcome),
797    /// The swing is under the group's evidence floor, where a phase is noise.
798    BelowFloor {
799        /// Peak-to-peak swing measured (metres).
800        amplitude_m: f64,
801        /// The floor it is under (metres).
802        floor_m: f64,
803    },
804    /// The measured stride anchor, in cycle fraction `[0, 1)`.
805    Measured(f64),
806}
807
808impl GaitMemberPhase {
809    /// The anchor a member contributes to its group's spread, if any.
810    pub fn anchor(self) -> Option<f64> {
811        match self {
812            GaitMemberPhase::Measured(phase) => Some(phase),
813            _ => None,
814        }
815    }
816}
817
818/// Classify one declared gait-group member.
819///
820/// `evidence` is what [`gait_phase_evidence`] returned for the member's clip,
821/// or `None` where the document holds no such clip (`present` is `false`) or
822/// the clip carries no readable foot cycle.
823pub fn gait_member_phase(
824    roles: &ResolvedRoles,
825    present: bool,
826    evidence: Option<&GaitPhaseEvidence>,
827    min_lr_amplitude_m: f64,
828) -> GaitMemberPhase {
829    if !present {
830        return GaitMemberPhase::Absent;
831    }
832    if !gait_roles_resolved(roles) {
833        return GaitMemberPhase::RolesUnresolved;
834    }
835    let Some(evidence) = evidence else {
836        return GaitMemberPhase::NoFootCycle;
837    };
838    match evidence.outcome {
839        GaitPhaseOutcome::Measured(_) if evidence.lr_amplitude_m < min_lr_amplitude_m => {
840            GaitMemberPhase::BelowFloor {
841                amplitude_m: evidence.lr_amplitude_m,
842                floor_m: min_lr_amplitude_m,
843            }
844        }
845        GaitPhaseOutcome::Measured(phase) => GaitMemberPhase::Measured(phase),
846        outcome => GaitMemberPhase::NoAnchor(outcome),
847    }
848}
849
850/// Normalized cycle position `[0,1)` of the minimum of the signal's
851/// first Fourier harmonic. Robust to plateaus and per-frame noise: the
852/// minimum of `A·cos(2π·t/N − φ)` sits at `t/N = (φ/2π + 0.5) mod 1`.
853pub fn fundamental_trough_phase(signal: &[f64]) -> Option<f64> {
854    let n = signal.len();
855    if n < 2 || signal.iter().any(|value| !value.is_finite()) {
856        return None;
857    }
858    let mut re = 0.0f64;
859    let mut im = 0.0f64;
860    for (k, y) in signal.iter().enumerate() {
861        let angle = std::f64::consts::TAU * k as f64 / n as f64;
862        re += y * angle.cos();
863        im += y * angle.sin();
864    }
865    let phi = im.atan2(re);
866    let phase = (phi / std::f64::consts::TAU + 0.5).rem_euclid(1.0);
867    phase.is_finite().then_some(phase)
868}
869
870/// Horizontal (XZ-plane) root displacement over the clip, divided by
871/// duration. Uses the Root role whenever it resolves and falls back to Hips
872/// only when Root is unresolved (clips without a dedicated root bone carry
873/// travel on the hips). Returns `None` rather than falling back when the
874/// selected role index is outside `grid`, the grid is too short, duration is
875/// non-positive, or the derived speed is non-finite.
876pub fn root_motion_speed_mps(grid: &PoseGrid, roles: &ResolvedRoles) -> Option<f64> {
877    let bone = roles.get(Role::Root).or_else(|| roles.get(Role::Hips))?;
878    let frames = grid.frame_count();
879    if frames < 2 || bone >= grid.bone_count() {
880        return None;
881    }
882    let duration = *grid.times.last()? as f64;
883    if duration <= 0.0 {
884        return None;
885    }
886    let a = grid.model_position(0, bone);
887    let b = grid.model_position(frames - 1, bone);
888    let dx = (b.x - a.x) as f64;
889    let dz = (b.z - a.z) as f64;
890    let speed = dx.hypot(dz) / duration;
891    speed.is_finite().then_some(speed)
892}
893
894/// Maximum angular deviation (degrees) of a rotation track from its
895/// first keyed rotation.
896pub fn rotation_range_deg(track: &Track) -> Option<f64> {
897    if track.property != Property::Rotation {
898        return None;
899    }
900    let first = track.key_quat(0)?;
901    if !first.is_finite() || first.length_squared() == 0.0 {
902        return None;
903    }
904    let first = first.normalize();
905    let mut max_deg = 0.0f64;
906    for k in 1..track.key_count() {
907        if let Some(q) = track.key_quat(k)
908            && q.is_finite()
909            && q.length_squared() > 0.0
910        {
911            let deg = first.angle_between(q.normalize()).to_degrees() as f64;
912            if deg.is_finite() {
913                max_deg = max_deg.max(deg);
914            }
915        }
916    }
917    Some(max_deg)
918}
919
920/// A cycle position folded into `[0, 1)`.
921///
922/// `rem_euclid` rounds a hair-negative input's wrap up to exactly one, and
923/// the cycle position of a full turn is zero.
924pub fn wrap_unit_phase(phase: f64) -> f64 {
925    let wrapped = phase.rem_euclid(1.0);
926    if wrapped < 1.0 { wrapped } else { 0.0 }
927}
928
929/// Circular distance between two normalized phases, in cycle fraction
930/// `[0, 0.5]`. Phases live on a ring, so the distance is the shorter of the
931/// two arcs between them.
932pub fn circular_phase_distance(phase: f64, other: f64) -> f64 {
933    let distance = (phase - other).abs() % 1.0;
934    distance.min(1.0 - distance)
935}
936
937/// The circular mean of a set of normalized phases and the maximum deviation
938/// from it, in one pass; `None` for an empty set, which has no mean
939/// direction.
940///
941/// The mean is a cycle position in `[0, 1)` and the deviation a cycle
942/// fraction in `[0, 0.5]`. A caller drawing a tolerance around a group of
943/// phases needs the centre as well as the spread, and taking them separately
944/// computed the resultant twice.
945///
946/// A set whose vectors cancel keeps whatever direction `atan2` derives from
947/// the residual sums — the same direction the spread has always been
948/// measured from.
949pub fn circular_phase_center_spread(phases: &[f64]) -> Option<(f64, f64)> {
950    use std::f64::consts::{PI, TAU};
951    if phases.is_empty() {
952        return None;
953    }
954    let (mut sin_sum, mut cos_sum) = (0.0f64, 0.0f64);
955    for p in phases {
956        sin_sum += (p * TAU).sin();
957        cos_sum += (p * TAU).cos();
958    }
959    let mean = sin_sum.atan2(cos_sum);
960    let mut max_dev = 0.0f64;
961    for p in phases {
962        let mut d = (p * TAU - mean).abs() % TAU;
963        if d > PI {
964            d = TAU - d;
965        }
966        max_dev = max_dev.max(d / TAU);
967    }
968    Some((wrap_unit_phase(mean / TAU), max_dev))
969}
970
971/// Maximum circular distance (in cycle fraction, `[0, 0.5]`) of a set of
972/// normalized phases from their circular mean. Phases live on a ring, so
973/// a naive max−min would over-report a cluster straddling the 0/1 wrap.
974pub fn circular_phase_spread(phases: &[f64]) -> f64 {
975    circular_phase_center_spread(phases).map_or(0.0, |(_, spread)| spread)
976}
977
978/// The metric sampling grid for a clip: uniform, resolution = max key
979/// count (mirroring how the runtime loops a clip over `[0, duration]`,
980/// wrapping duration→0 at render times unaligned with authored keys).
981/// `None` for clips too short to carry a cycle (< 3 keys), matching the
982/// reference implementation.
983pub fn metric_frame_count(clip: &Clip) -> Option<usize> {
984    let n = crate::sample::default_frame_count(clip);
985    if clip.duration_s <= 0.0 || n < 3 {
986        None
987    } else {
988        Some(n)
989    }
990}
991
992#[cfg(test)]
993mod tests {
994    use super::*;
995    use crate::check::CheckCtx;
996    use crate::config::Config;
997    use crate::measure::{RootTrajectorySourceRole, measure_document};
998    use crate::model::{
999        Bone, Clip, Document, Interpolation, Property, Skeleton, Track, TrackValues, Transform,
1000    };
1001    use crate::profile::{ResolvedRoles, Role};
1002    use glam::{EulerRot, Mat3, Quat, Vec3};
1003    use std::rc::Rc;
1004
1005    #[test]
1006    fn gait_phase_outcome_retains_the_defensive_derivation_failure_state() {
1007        assert_eq!(
1008            GaitPhaseOutcome::classify(None, 0.01, true),
1009            GaitPhaseOutcome::Unavailable
1010        );
1011        assert_eq!(
1012            GaitPhaseOutcome::classify(Some(0.25), 0.01, true),
1013            GaitPhaseOutcome::Measured(0.25)
1014        );
1015        assert_eq!(
1016            GaitPhaseOutcome::classify(None, 0.0, true),
1017            GaitPhaseOutcome::NoFootHeightSwing
1018        );
1019        assert_eq!(
1020            GaitPhaseOutcome::classify(None, 0.0, false),
1021            GaitPhaseOutcome::MissingBilateralFootRoles
1022        );
1023    }
1024
1025    fn document_with_metric_clip() -> Document {
1026        Document {
1027            skeleton: Skeleton {
1028                bones: vec![Bone {
1029                    name: "root".into(),
1030                    parent: None,
1031                    rest: Transform::IDENTITY,
1032                    inverse_bind: None,
1033                }],
1034            },
1035            clips: vec![Clip {
1036                name: "walk".into(),
1037                duration_s: 1.0,
1038                tracks: vec![Track {
1039                    bone: 0,
1040                    property: Property::Rotation,
1041                    interpolation: Interpolation::Linear,
1042                    times: vec![0.0, 0.5, 1.0],
1043                    values: TrackValues::Quats(vec![
1044                        Quat::IDENTITY,
1045                        Quat::from_rotation_y(0.1),
1046                        Quat::from_rotation_y(0.2),
1047                    ]),
1048                }],
1049            }],
1050            ..Document::default()
1051        }
1052    }
1053
1054    fn document_with_grid_inputs(duration_s: f64, times: Vec<f32>) -> Document {
1055        let values = vec![Quat::IDENTITY; times.len()];
1056        Document {
1057            skeleton: Skeleton {
1058                bones: vec![Bone {
1059                    name: "root".into(),
1060                    parent: None,
1061                    rest: Transform::IDENTITY,
1062                    inverse_bind: None,
1063                }],
1064            },
1065            clips: vec![Clip {
1066                name: "probe".into(),
1067                duration_s,
1068                tracks: vec![Track {
1069                    bone: 0,
1070                    property: Property::Rotation,
1071                    interpolation: Interpolation::Linear,
1072                    times,
1073                    values: TrackValues::Quats(values),
1074                }],
1075            }],
1076            ..Document::default()
1077        }
1078    }
1079
1080    #[test]
1081    fn metric_grids_are_shared_by_checks_and_measurements() {
1082        let doc = document_with_metric_clip();
1083        let roles = ResolvedRoles::default();
1084        let config = Config::default();
1085        let grids = MetricGrids::new(&doc);
1086
1087        let ctx = CheckCtx::new(&grids, &roles, &config);
1088        let from_ctx = ctx.grid(0).expect("metric grid");
1089        let from_owner = grids.grid(0).expect("same metric grid");
1090        assert!(Rc::ptr_eq(&from_ctx, &from_owner));
1091
1092        let measurements = measure_document(&grids, &roles, &config);
1093        assert!(measurements.contains_key("walk"));
1094        let fresh_grids = MetricGrids::new(&doc);
1095        assert_eq!(
1096            serde_json::to_value(&measurements).expect("shared measurements serialize"),
1097            serde_json::to_value(measure_document(&fresh_grids, &roles, &config))
1098                .expect("plain measurements serialize")
1099        );
1100    }
1101
1102    #[test]
1103    fn grid_returns_none_for_each_documented_invalid_request() {
1104        let valid = document_with_grid_inputs(1.0, vec![0.0, 0.5, 1.0]);
1105        let valid_grids = MetricGrids::new(&valid);
1106        assert!(valid_grids.grid(0).is_some());
1107        for clip_index in [1, 2, usize::MAX] {
1108            assert!(valid_grids.grid(clip_index).is_none());
1109        }
1110
1111        for duration_s in [0.0, -1.0] {
1112            let non_positive = document_with_grid_inputs(duration_s, vec![0.0, 0.5, 1.0]);
1113            assert!(MetricGrids::new(&non_positive).grid(0).is_none());
1114        }
1115
1116        for times in [vec![], vec![0.0], vec![0.0, 1.0]] {
1117            let too_few_keys = document_with_grid_inputs(1.0, times);
1118            assert!(MetricGrids::new(&too_few_keys).grid(0).is_none());
1119        }
1120    }
1121
1122    #[test]
1123    fn grid_uses_longest_track_for_resolution() {
1124        // The first track is too short by itself; the later translation
1125        // track selects the grid's three-frame resolution.
1126        let mut doc = document_with_grid_inputs(1.0, vec![0.0, 1.0]);
1127        doc.clips[0].tracks.push(Track {
1128            bone: 0,
1129            property: Property::Translation,
1130            interpolation: Interpolation::Linear,
1131            times: vec![0.0, 0.5, 1.0],
1132            values: TrackValues::Vec3s(vec![Vec3::ZERO, Vec3::X, 2.0 * Vec3::X]),
1133        });
1134
1135        let grid = MetricGrids::new(&doc)
1136            .grid(0)
1137            .expect("later longest track supplies a metric grid");
1138        assert_eq!(grid.frame_count(), 3);
1139    }
1140
1141    #[test]
1142    fn unrelated_dense_track_changes_shared_sampled_trajectory_not_the_analytic_curve() {
1143        let skeleton = Skeleton {
1144            bones: vec![
1145                Bone {
1146                    name: "root".into(),
1147                    parent: None,
1148                    rest: Transform::IDENTITY,
1149                    inverse_bind: None,
1150                },
1151                Bone {
1152                    name: "unrelated".into(),
1153                    parent: None,
1154                    rest: Transform::IDENTITY,
1155                    inverse_bind: None,
1156                },
1157            ],
1158        };
1159        let target_tracks = vec![
1160            Track {
1161                bone: 0,
1162                property: Property::Translation,
1163                interpolation: Interpolation::Linear,
1164                times: vec![0.0, 0.25, 1.0],
1165                values: TrackValues::Vec3s(vec![Vec3::ZERO, Vec3::new(1.0, 2.0, 0.0), Vec3::ZERO]),
1166            },
1167            Track {
1168                bone: 0,
1169                property: Property::Rotation,
1170                interpolation: Interpolation::Linear,
1171                times: vec![0.0, 0.25, 1.0],
1172                values: TrackValues::Quats(
1173                    [0.0_f32, 170.0, 340.0]
1174                        .map(|degrees| Quat::from_rotation_y(degrees.to_radians()))
1175                        .to_vec(),
1176                ),
1177            },
1178        ];
1179        let document = Document {
1180            skeleton: skeleton.clone(),
1181            clips: vec![Clip {
1182                name: "aliased".into(),
1183                duration_s: 1.0,
1184                tracks: target_tracks.clone(),
1185            }],
1186            ..Document::default()
1187        };
1188        let coarse_grid = MetricGrids::new(&document)
1189            .grid(0)
1190            .expect("three-sample grid");
1191        assert_eq!(coarse_grid.frame_count(), 3);
1192        let coarse = root_trajectory_metrics(&coarse_grid, 0).expect("coarse trajectory");
1193        let coarse_translation = coarse.translation.expect("coarse translation");
1194        let coarse_yaw = coarse.yaw.expect("coarse yaw");
1195
1196        let analytic_peak_y = 2.0;
1197        assert!((coarse_translation.horizontal_travel_m - 4.0 / 3.0).abs() < 1.0e-5);
1198        assert!((coarse_translation.vertical_max_displacement_m - 4.0 / 3.0).abs() < 1.0e-5);
1199        assert!(coarse_translation.vertical_max_displacement_m < analytic_peak_y);
1200        assert!((coarse_yaw.net_yaw_deg - -20.0).abs() < 1.0e-4);
1201        assert!((coarse_yaw.unwrapped_yaw_deg - -20.0).abs() < 1.0e-4);
1202        assert!((coarse_yaw.yaw_travel_deg - 740.0 / 3.0).abs() < 1.0e-3);
1203
1204        let mut dense_document = Document {
1205            skeleton,
1206            clips: vec![Clip {
1207                name: "aliased".into(),
1208                duration_s: 1.0,
1209                tracks: target_tracks,
1210            }],
1211            ..Document::default()
1212        };
1213        dense_document.clips[0].tracks.push(Track {
1214            bone: 1,
1215            property: Property::Scale,
1216            interpolation: Interpolation::Linear,
1217            times: vec![0.0, 0.25, 0.5, 0.75, 1.0],
1218            values: TrackValues::Vec3s(vec![Vec3::ONE; 5]),
1219        });
1220        let dense_grid = MetricGrids::new(&dense_document)
1221            .grid(0)
1222            .expect("unrelated track selects five-sample grid");
1223        assert_eq!(dense_grid.frame_count(), 5);
1224        let dense = root_trajectory_metrics(&dense_grid, 0).expect("dense trajectory");
1225        let dense_translation = dense.translation.expect("dense translation");
1226        let dense_yaw = dense.yaw.expect("dense yaw");
1227
1228        assert!((dense_translation.horizontal_travel_m - 2.0).abs() < 1.0e-5);
1229        assert_eq!(
1230            dense_translation.vertical_max_displacement_m,
1231            analytic_peak_y
1232        );
1233        assert!((dense_yaw.net_yaw_deg - -20.0).abs() < 1.0e-4);
1234        assert!((dense_yaw.unwrapped_yaw_deg - 340.0).abs() < 1.0e-3);
1235        assert!((dense_yaw.yaw_travel_deg - 340.0).abs() < 1.0e-3);
1236
1237        let coarse_roles =
1238            ResolvedRoles::from_names(&document.skeleton, [(Role::Root, "root".into())]);
1239        let coarse_measurements = measure_document(
1240            &MetricGrids::new(&document),
1241            &coarse_roles,
1242            &Config::default(),
1243        );
1244        let coarse_published = coarse_measurements["aliased"]
1245            .root_trajectory
1246            .as_ref()
1247            .expect("public coarse Root trajectory");
1248        assert_eq!(coarse_published.source_role, RootTrajectorySourceRole::Root);
1249        let coarse_published_translation = coarse_published.translation.unwrap();
1250        let coarse_published_yaw = coarse_published.yaw.unwrap();
1251        assert!((coarse_published_translation.horizontal_travel_m - 4.0 / 3.0).abs() < 1.0e-5);
1252        assert!(
1253            (coarse_published_translation.vertical_max_displacement_m - 4.0 / 3.0).abs() < 1.0e-5
1254        );
1255        assert!((coarse_published_yaw.unwrapped_yaw_deg - -20.0).abs() < 1.0e-4);
1256        assert!((coarse_published_yaw.yaw_travel_deg - 740.0 / 3.0).abs() < 1.0e-3);
1257
1258        let dense_roles =
1259            ResolvedRoles::from_names(&dense_document.skeleton, [(Role::Root, "root".into())]);
1260        let dense_measurements = measure_document(
1261            &MetricGrids::new(&dense_document),
1262            &dense_roles,
1263            &Config::default(),
1264        );
1265        let dense_published = dense_measurements["aliased"]
1266            .root_trajectory
1267            .as_ref()
1268            .expect("public dense Root trajectory");
1269        assert_eq!(dense_published.source_role, RootTrajectorySourceRole::Root);
1270        let dense_published_translation = dense_published.translation.unwrap();
1271        let dense_published_yaw = dense_published.yaw.unwrap();
1272        assert!((dense_published_translation.horizontal_travel_m - 2.0).abs() < 1.0e-5);
1273        assert_eq!(
1274            dense_published_translation.vertical_max_displacement_m,
1275            analytic_peak_y
1276        );
1277        assert!((dense_published_yaw.unwrapped_yaw_deg - 340.0).abs() < 1.0e-3);
1278        assert!((dense_published_yaw.yaw_travel_deg - 340.0).abs() < 1.0e-3);
1279        assert_ne!(
1280            coarse_published_translation.horizontal_travel_m,
1281            dense_published_translation.horizontal_travel_m
1282        );
1283        assert_ne!(
1284            coarse_published_translation.vertical_max_displacement_m,
1285            dense_published_translation.vertical_max_displacement_m
1286        );
1287        assert_ne!(
1288            coarse_published_yaw.unwrapped_yaw_deg,
1289            dense_published_yaw.unwrapped_yaw_deg
1290        );
1291    }
1292
1293    #[test]
1294    fn foot_metrics_reject_finite_positions_whose_relative_subtraction_overflows() {
1295        let mut doc = document_with_metric_clip();
1296        doc.skeleton.bones = vec![
1297            Bone {
1298                name: "hips".into(),
1299                parent: None,
1300                rest: Transform {
1301                    translation: Vec3::splat(-f32::MAX),
1302                    ..Transform::IDENTITY
1303                },
1304                inverse_bind: None,
1305            },
1306            Bone {
1307                name: "left".into(),
1308                parent: None,
1309                rest: Transform {
1310                    translation: Vec3::splat(f32::MAX),
1311                    ..Transform::IDENTITY
1312                },
1313                inverse_bind: None,
1314            },
1315        ];
1316        doc.clips[0].tracks[0].bone = 0;
1317        let roles = ResolvedRoles::from_names(
1318            &doc.skeleton,
1319            [
1320                (Role::Hips, "hips".to_string()),
1321                (Role::LeftFoot, "left".to_string()),
1322            ],
1323        );
1324        let grid = MetricGrids::new(&doc).grid(0).expect("metric grid");
1325
1326        assert!(grid.model_position(0, 0).is_finite());
1327        assert!(grid.model_position(0, 1).is_finite());
1328        assert!(foot_cycle_metrics(&grid, &roles, MIN_STRIDE_STEP_M).is_none());
1329    }
1330
1331    /// A real stride can still leave `loop_seam_ratio` `None`: the seam
1332    /// (frame `frames - 1` vs frame `0`) is a single point-to-point
1333    /// distance, unconstrained by either neighbour step, so it can sit at
1334    /// a per-axis delta near `f32::MAX` while both neighbour steps stay at
1335    /// the smallest representable positive `f32` value (comfortably over
1336    /// the configured floor). `f32` squares that delta while computing the
1337    /// distance, overflowing to infinity even though every input position
1338    /// was finite — the ratio then divides out non-finite, not "no
1339    /// subject". This is the only known route to
1340    /// [`FootCycleMetrics::loop_seam_ratio`]'s "real stride but
1341    /// underivable" `None`, and it requires magnitudes far outside any
1342    /// real animation.
1343    #[test]
1344    fn foot_metrics_real_stride_with_seam_beyond_f32_squaring_range_has_no_ratio() {
1345        let mut doc = document_with_metric_clip();
1346        doc.skeleton.bones = vec![
1347            Bone {
1348                name: "hips".into(),
1349                parent: None,
1350                rest: Transform::IDENTITY,
1351                inverse_bind: None,
1352            },
1353            Bone {
1354                name: "left".into(),
1355                parent: Some(0),
1356                rest: Transform::IDENTITY,
1357                inverse_bind: None,
1358            },
1359        ];
1360        doc.clips[0].tracks = vec![Track {
1361            bone: 1,
1362            property: Property::Translation,
1363            interpolation: Interpolation::Linear,
1364            times: vec![0.0, 0.25, 0.5, 1.0],
1365            values: TrackValues::Vec3s(vec![
1366                Vec3::ZERO,
1367                Vec3::new(f32::MIN_POSITIVE, 0.0, 0.0),
1368                Vec3::new(f32::MAX - f32::MIN_POSITIVE, 0.0, 0.0),
1369                Vec3::new(f32::MAX, 0.0, 0.0),
1370            ]),
1371        }];
1372        let roles = ResolvedRoles::from_names(
1373            &doc.skeleton,
1374            [
1375                (Role::Hips, "hips".to_string()),
1376                (Role::LeftFoot, "left".to_string()),
1377            ],
1378        );
1379        let grid = MetricGrids::new(&doc).grid(0).expect("metric grid");
1380
1381        let metrics = foot_cycle_metrics(&grid, &roles, f64::from(f32::MIN_POSITIVE))
1382            .expect("hips + one foot role with enough frames yields metrics");
1383
1384        assert!(
1385            metrics.has_real_stride,
1386            "the neighbour step met the (tiny) configured floor"
1387        );
1388        assert_eq!(
1389            metrics.loop_seam_ratio, None,
1390            "the seam distance overflowed f32 squaring to infinity, so the \
1391             ratio is non-finite despite a real stride"
1392        );
1393    }
1394
1395    fn trajectory_grid_with_rotations(positions: Vec<Vec3>, rotations: Vec<Quat>) -> PoseGrid {
1396        assert_eq!(positions.len(), rotations.len());
1397        assert!(positions.len() >= 2);
1398        let last = positions.len() - 1;
1399        let times = (0..=last)
1400            .map(|index| index as f32 / last as f32)
1401            .collect::<Vec<_>>();
1402        let skeleton = Skeleton {
1403            bones: vec![Bone {
1404                name: "root".into(),
1405                parent: None,
1406                rest: Transform::IDENTITY,
1407                inverse_bind: None,
1408            }],
1409        };
1410        let clip = Clip {
1411            name: "trajectory".into(),
1412            duration_s: 1.0,
1413            tracks: vec![
1414                Track {
1415                    bone: 0,
1416                    property: Property::Translation,
1417                    interpolation: Interpolation::Linear,
1418                    times: times.clone(),
1419                    values: TrackValues::Vec3s(positions),
1420                },
1421                Track {
1422                    bone: 0,
1423                    property: Property::Rotation,
1424                    interpolation: Interpolation::Linear,
1425                    times,
1426                    values: TrackValues::Quats(rotations),
1427                },
1428            ],
1429        };
1430        sample_clip(&skeleton, &clip, last + 1)
1431    }
1432
1433    fn trajectory_grid(positions: Vec<Vec3>, yaw_degrees: Vec<f32>) -> PoseGrid {
1434        trajectory_grid_with_rotations(
1435            positions,
1436            yaw_degrees
1437                .into_iter()
1438                .map(|degrees| Quat::from_rotation_y(degrees.to_radians()))
1439                .collect(),
1440        )
1441    }
1442
1443    fn assert_stationary_translation(trajectory: &RootTrajectoryMetrics) {
1444        let translation = trajectory
1445            .translation
1446            .expect("finite stationary translation");
1447        assert_eq!(translation.horizontal_displacement_x_m, 0.0);
1448        assert_eq!(translation.horizontal_displacement_z_m, 0.0);
1449        assert_eq!(translation.horizontal_travel_m, 0.0);
1450        assert_eq!(translation.vertical_displacement_m, 0.0);
1451        assert_eq!(translation.vertical_min_displacement_m, 0.0);
1452        assert_eq!(translation.vertical_max_displacement_m, 0.0);
1453    }
1454
1455    #[test]
1456    fn root_trajectory_retains_direction_travel_and_vertical_extrema() {
1457        let grid = trajectory_grid(
1458            vec![
1459                Vec3::ZERO,
1460                Vec3::new(1.0, 2.0, 0.0),
1461                Vec3::new(1.0, 1.0, -1.0),
1462            ],
1463            vec![0.0, 45.0, 90.0],
1464        );
1465        let trajectory = root_trajectory_metrics(&grid, 0).expect("valid selected bone");
1466        let translation = trajectory.translation.expect("finite translation");
1467
1468        assert_eq!(translation.horizontal_displacement_x_m, 1.0);
1469        assert_eq!(translation.horizontal_displacement_z_m, -1.0);
1470        assert_eq!(translation.horizontal_travel_m, 2.0);
1471        assert!(
1472            translation.horizontal_travel_m
1473                > translation
1474                    .horizontal_displacement_x_m
1475                    .hypot(translation.horizontal_displacement_z_m)
1476        );
1477        assert_eq!(translation.vertical_displacement_m, 1.0);
1478        assert_eq!(translation.vertical_min_displacement_m, 0.0);
1479        assert_eq!(translation.vertical_max_displacement_m, 2.0);
1480        let yaw = trajectory.yaw.expect("finite yaw");
1481        assert!((yaw.net_yaw_deg - 90.0).abs() < 1.0e-4);
1482        assert!((yaw.unwrapped_yaw_deg - 90.0).abs() < 1.0e-4);
1483
1484        let out_and_back =
1485            trajectory_grid(vec![Vec3::ZERO, Vec3::X, Vec3::ZERO], vec![0.0, 0.0, 0.0]);
1486        let translation = root_trajectory_metrics(&out_and_back, 0)
1487            .unwrap()
1488            .translation
1489            .unwrap();
1490        assert_eq!(translation.horizontal_displacement_x_m, 0.0);
1491        assert_eq!(translation.horizontal_travel_m, 2.0);
1492
1493        for (positions, expected_x, expected_z, expected_travel) in [
1494            (vec![Vec3::ZERO, -Vec3::Z, -Vec3::Z * 2.0], 0.0, -2.0, 2.0),
1495            (vec![Vec3::ZERO, Vec3::X * 0.5, Vec3::X], 1.0, 0.0, 1.0),
1496        ] {
1497            let grid = trajectory_grid(positions, vec![0.0; 3]);
1498            let translation = root_trajectory_metrics(&grid, 0)
1499                .unwrap()
1500                .translation
1501                .unwrap();
1502            assert_eq!(translation.horizontal_displacement_x_m, expected_x);
1503            assert_eq!(translation.horizontal_displacement_z_m, expected_z);
1504            assert_eq!(translation.horizontal_travel_m, expected_travel);
1505        }
1506
1507        for (positions, expected_net, expected_min, expected_max) in [
1508            (vec![Vec3::ZERO, Vec3::Y, Vec3::Y * 2.0], 2.0, 0.0, 2.0),
1509            (vec![Vec3::ZERO, -Vec3::Y, -Vec3::Y * 2.0], -2.0, -2.0, 0.0),
1510            (vec![Vec3::ZERO, Vec3::Y * 2.0, Vec3::ZERO], 0.0, 0.0, 2.0),
1511        ] {
1512            let grid = trajectory_grid(positions, vec![0.0; 3]);
1513            let translation = root_trajectory_metrics(&grid, 0)
1514                .unwrap()
1515                .translation
1516                .unwrap();
1517            assert_eq!(translation.vertical_displacement_m, expected_net);
1518            assert_eq!(translation.vertical_min_displacement_m, expected_min);
1519            assert_eq!(translation.vertical_max_displacement_m, expected_max);
1520        }
1521    }
1522
1523    #[test]
1524    fn root_yaw_retains_signed_half_and_full_turns_and_reversing_travel() {
1525        assert_eq!(canonical_net_yaw_deg(179.99995), 180.0);
1526        assert_eq!(canonical_net_yaw_deg(-179.99995), -180.0);
1527        let cases = [
1528            (vec![0.0, 45.0, 90.0], 90.0, 90.0, 90.0),
1529            (vec![0.0, -45.0, -90.0], -90.0, -90.0, 90.0),
1530            (vec![0.0, 90.0, 180.0], 180.0, 180.0, 180.0),
1531            (vec![0.0, -90.0, -180.0], -180.0, -180.0, 180.0),
1532            (
1533                vec![0.0, 60.0, 120.0, 179.99995],
1534                180.0,
1535                179.99995,
1536                179.99995,
1537            ),
1538            (
1539                vec![0.0, -60.0, -120.0, -179.99995],
1540                -180.0,
1541                -179.99995,
1542                179.99995,
1543            ),
1544            (vec![0.0, 90.0, 180.0, 270.0, 360.0], 0.0, 360.0, 360.0),
1545            (vec![0.0, -90.0, -180.0, -270.0, -360.0], 0.0, -360.0, 360.0),
1546            (vec![0.0, 90.0, 0.0], 0.0, 0.0, 180.0),
1547        ];
1548
1549        for (angles, expected_net, expected_unwrapped, expected_travel) in cases {
1550            let grid = trajectory_grid(vec![Vec3::ZERO; angles.len()], angles);
1551            let trajectory = root_trajectory_metrics(&grid, 0).unwrap();
1552            assert_stationary_translation(&trajectory);
1553            let yaw = trajectory.yaw.expect("yaw with sub-half-turn steps");
1554            assert_eq!(yaw.heading_axis, RootYawHeadingAxis::PositiveZ);
1555            if expected_net == 180.0 || expected_net == -180.0 {
1556                assert_eq!(yaw.net_yaw_deg, expected_net, "{yaw:?}");
1557            } else {
1558                assert!((yaw.net_yaw_deg - expected_net).abs() < 1.0e-4, "{yaw:?}");
1559            }
1560            assert!(
1561                (yaw.unwrapped_yaw_deg - expected_unwrapped).abs() < 1.0e-4,
1562                "{yaw:?}"
1563            );
1564            assert!(
1565                (yaw.yaw_travel_deg - expected_travel).abs() < 1.0e-4,
1566                "{yaw:?}"
1567            );
1568        }
1569    }
1570
1571    #[test]
1572    fn root_trajectory_translation_and_yaw_fail_independently() {
1573        let bad_position = trajectory_grid(
1574            vec![Vec3::ZERO, Vec3::new(f32::NAN, 0.0, 0.0), Vec3::ZERO],
1575            vec![0.0, 45.0, 90.0],
1576        );
1577        let measured = root_trajectory_metrics(&bad_position, 0).unwrap();
1578        assert!(measured.translation.is_none());
1579        assert!(measured.yaw.is_some());
1580
1581        let bad_rotation = trajectory_grid_with_rotations(
1582            vec![Vec3::ZERO, Vec3::X, Vec3::X * 2.0],
1583            vec![
1584                Quat::IDENTITY,
1585                Quat::from_xyzw(f32::NAN, 0.0, 0.0, 1.0),
1586                Quat::IDENTITY,
1587            ],
1588        );
1589        let measured = root_trajectory_metrics(&bad_rotation, 0).unwrap();
1590        let translation = measured
1591            .translation
1592            .expect("translation remains measurable");
1593        assert_eq!(translation.horizontal_displacement_x_m, 2.0);
1594        assert_eq!(translation.horizontal_displacement_z_m, 0.0);
1595        assert_eq!(translation.horizontal_travel_m, 2.0);
1596        assert_eq!(translation.vertical_displacement_m, 0.0);
1597        assert_eq!(translation.vertical_min_displacement_m, 0.0);
1598        assert_eq!(translation.vertical_max_displacement_m, 0.0);
1599        assert!(measured.yaw.is_none());
1600
1601        let zero_rotation = trajectory_grid_with_rotations(
1602            vec![Vec3::ZERO, Vec3::X, Vec3::X * 2.0],
1603            vec![
1604                Quat::IDENTITY,
1605                Quat::from_xyzw(0.0, 0.0, 0.0, 0.0),
1606                Quat::IDENTITY,
1607            ],
1608        );
1609        let measured = root_trajectory_metrics(&zero_rotation, 0).unwrap();
1610        let translation = measured
1611            .translation
1612            .expect("zero quaternion only makes rotation decomposition unavailable");
1613        assert_eq!(translation.horizontal_displacement_x_m, 2.0);
1614        assert_eq!(translation.horizontal_displacement_z_m, 0.0);
1615        assert_eq!(translation.horizontal_travel_m, 2.0);
1616        assert_eq!(translation.vertical_displacement_m, 0.0);
1617        assert_eq!(translation.vertical_min_displacement_m, 0.0);
1618        assert_eq!(translation.vertical_max_displacement_m, 0.0);
1619        assert!(measured.yaw.is_none());
1620    }
1621
1622    #[test]
1623    fn root_trajectory_extrema_widen_finite_samples_and_reject_model_overflow() {
1624        let finite_extremes = trajectory_grid(
1625            vec![-Vec3::Y * f32::MAX, Vec3::Y * f32::MAX],
1626            vec![0.0, 0.0],
1627        );
1628        let translation = root_trajectory_metrics(&finite_extremes, 0)
1629            .expect("selected bone exists")
1630            .translation
1631            .expect("finite binary32 samples have finite widened extrema");
1632        let full_binary32_span = 2.0 * f64::from(f32::MAX);
1633        assert!(full_binary32_span.is_finite());
1634        assert_eq!(translation.vertical_displacement_m, full_binary32_span);
1635        assert_eq!(translation.vertical_min_displacement_m, 0.0);
1636        assert_eq!(translation.vertical_max_displacement_m, full_binary32_span);
1637
1638        let skeleton = Skeleton {
1639            bones: vec![
1640                Bone {
1641                    name: "ancestor".into(),
1642                    parent: None,
1643                    rest: Transform {
1644                        translation: Vec3::Y * f32::MAX,
1645                        ..Transform::IDENTITY
1646                    },
1647                    inverse_bind: None,
1648                },
1649                Bone {
1650                    name: "root".into(),
1651                    parent: Some(0),
1652                    rest: Transform {
1653                        translation: Vec3::Y * f32::MAX,
1654                        ..Transform::IDENTITY
1655                    },
1656                    inverse_bind: None,
1657                },
1658            ],
1659        };
1660        let clip = Clip {
1661            name: "overflow".into(),
1662            duration_s: 1.0,
1663            tracks: vec![Track {
1664                bone: 1,
1665                property: Property::Rotation,
1666                interpolation: Interpolation::Linear,
1667                times: vec![0.0, 1.0],
1668                values: TrackValues::Quats(vec![Quat::IDENTITY; 2]),
1669            }],
1670        };
1671        let grid = sample_clip(&skeleton, &clip, 2);
1672
1673        assert!(
1674            skeleton
1675                .bones
1676                .iter()
1677                .all(|bone| bone.rest.translation.is_finite()
1678                    && bone.rest.rotation.is_finite()
1679                    && bone.rest.scale.is_finite())
1680        );
1681        assert!(
1682            !grid.model_position(0, 1).is_finite(),
1683            "finite local translations overflow while composing model space"
1684        );
1685        let trajectory = root_trajectory_metrics(&grid, 1).expect("selected bone exists");
1686        assert_eq!(trajectory.translation, None);
1687        let yaw = trajectory
1688            .yaw
1689            .expect("translation overflow does not erase finite yaw");
1690        assert_eq!(yaw.net_yaw_deg, 0.0);
1691        assert_eq!(yaw.unwrapped_yaw_deg, 0.0);
1692        assert_eq!(yaw.yaw_travel_deg, 0.0);
1693    }
1694
1695    #[test]
1696    fn root_yaw_publishes_fixed_basis_and_refuses_ambiguous_steps() {
1697        let in_place_grid =
1698            trajectory_grid(vec![Vec3::new(3.0, 2.0, -1.0); 3], vec![0.0, 45.0, 90.0]);
1699        let in_place = root_trajectory_metrics(&in_place_grid, 0).unwrap();
1700        assert_stationary_translation(&in_place);
1701
1702        let positive_x_rotation = Quat::from_rotation_x(std::f32::consts::FRAC_PI_6);
1703        let positive_x_length = {
1704            let (x, z) = horizontal_heading(
1705                positive_x_rotation.as_dquat(),
1706                RootYawHeadingAxis::PositiveX,
1707            );
1708            x.hypot(z)
1709        };
1710        for other_axis in [RootYawHeadingAxis::PositiveZ, RootYawHeadingAxis::PositiveY] {
1711            let (x, z) = horizontal_heading(positive_x_rotation.as_dquat(), other_axis);
1712            assert!(positive_x_length > x.hypot(z));
1713        }
1714        let positive_x_grid =
1715            trajectory_grid_with_rotations(vec![Vec3::ZERO; 3], vec![positive_x_rotation; 3]);
1716        let positive_x = root_trajectory_metrics(&positive_x_grid, 0).unwrap();
1717        assert_stationary_translation(&positive_x);
1718        assert_eq!(
1719            positive_x.yaw.unwrap().heading_axis,
1720            RootYawHeadingAxis::PositiveX
1721        );
1722
1723        // Analytic tilted basis, built directly from orthonormal columns rather
1724        // than Euler angles. Local +Y has world heading 30 degrees and a 0.5
1725        // vertical component; +X and +Z each have a smaller horizontal
1726        // projection, so +Y is the strict fixed witness.
1727        let sqrt_three = 3.0_f32.sqrt();
1728        let positive_y = Vec3::new(sqrt_three / 4.0, 0.5, 0.75);
1729        let horizontal_perpendicular = Vec3::new(sqrt_three / 2.0, 0.0, -0.5);
1730        let tilt_tangent = Vec3::new(-0.25, sqrt_three / 2.0, -sqrt_three / 4.0);
1731        let positive_x =
1732            (horizontal_perpendicular + tilt_tangent) * std::f32::consts::FRAC_1_SQRT_2;
1733        let positive_z = positive_x.cross(positive_y);
1734        let tilted = Quat::from_mat3(&Mat3::from_cols(positive_x, positive_y, positive_z));
1735        let (heading_x, heading_z) =
1736            horizontal_heading(tilted.as_dquat(), RootYawHeadingAxis::PositiveY);
1737        assert!((heading_x.atan2(heading_z).to_degrees() - 30.0).abs() < 1.0e-5);
1738        let rotations = [0.0_f32, 20.0, 40.0]
1739            .map(|degrees| Quat::from_axis_angle(Vec3::Y, degrees.to_radians()) * tilted)
1740            .to_vec();
1741        let basis_grid = trajectory_grid_with_rotations(vec![Vec3::ZERO; 3], rotations);
1742        let basis = root_trajectory_metrics(&basis_grid, 0).unwrap();
1743        assert_stationary_translation(&basis);
1744        let yaw = basis.yaw.unwrap();
1745        assert_eq!(yaw.heading_axis, RootYawHeadingAxis::PositiveY);
1746        assert!((yaw.net_yaw_deg - 40.0).abs() < 1.0e-4);
1747        assert!((yaw.unwrapped_yaw_deg - 40.0).abs() < 1.0e-4);
1748        assert!((yaw.yaw_travel_deg - 40.0).abs() < 1.0e-4);
1749
1750        // A local-Y twist changes a Y-first Euler decomposition but cannot
1751        // change the retained local +Y witness. This is the concrete case an
1752        // Euler-component shortcut gets wrong.
1753        let local_half_turn = Quat::from_xyzw(0.0, 1.0, 0.0, 0.0);
1754        let fixed_witness_rotations = vec![tilted, tilted * local_half_turn, tilted];
1755        let first_euler_y = fixed_witness_rotations[0].to_euler(EulerRot::YXZ).0;
1756        let twisted_euler_y = fixed_witness_rotations[1].to_euler(EulerRot::YXZ).0;
1757        assert!(
1758            (twisted_euler_y - first_euler_y).abs() > 1.0,
1759            "Y-first Euler shortcut must observe a misleading change"
1760        );
1761        let expected_heading = horizontal_heading(
1762            fixed_witness_rotations[0].as_dquat(),
1763            RootYawHeadingAxis::PositiveY,
1764        );
1765        for rotation in &fixed_witness_rotations[1..] {
1766            let heading = horizontal_heading(rotation.as_dquat(), RootYawHeadingAxis::PositiveY);
1767            assert!((heading.0 - expected_heading.0).abs() < 1.0e-7);
1768            assert!((heading.1 - expected_heading.1).abs() < 1.0e-7);
1769        }
1770
1771        let stationary_fixed_witness =
1772            trajectory_grid_with_rotations(vec![Vec3::ZERO; 3], fixed_witness_rotations.clone());
1773        let misleading_translation_fixed_witness = trajectory_grid_with_rotations(
1774            vec![Vec3::ZERO, Vec3::X * 100.0, Vec3::Z * 100.0],
1775            fixed_witness_rotations,
1776        );
1777        let stationary = root_trajectory_metrics(&stationary_fixed_witness, 0).unwrap();
1778        assert_stationary_translation(&stationary);
1779        let misleading = root_trajectory_metrics(&misleading_translation_fixed_witness, 0)
1780            .expect("translation content does not select yaw/up");
1781        assert!(misleading.translation.unwrap().horizontal_travel_m > 200.0);
1782        let stationary_yaw = stationary.yaw.unwrap();
1783        let misleading_yaw = misleading.yaw.unwrap();
1784        assert_eq!(stationary_yaw.heading_axis, RootYawHeadingAxis::PositiveY);
1785        assert_eq!(misleading_yaw, stationary_yaw);
1786        assert_eq!(stationary_yaw.net_yaw_deg, 0.0);
1787        assert_eq!(stationary_yaw.unwrapped_yaw_deg, 0.0);
1788        assert_eq!(stationary_yaw.yaw_travel_deg, 0.0);
1789
1790        let becomes_vertical = trajectory_grid_with_rotations(
1791            vec![Vec3::ZERO; 3],
1792            vec![
1793                Quat::IDENTITY,
1794                Quat::from_rotation_x(std::f32::consts::FRAC_PI_4),
1795                Quat::from_rotation_x(std::f32::consts::FRAC_PI_2),
1796            ],
1797        );
1798        assert!(
1799            root_trajectory_metrics(&becomes_vertical, 0)
1800                .unwrap()
1801                .yaw
1802                .is_none()
1803        );
1804
1805        let half_turn_step = trajectory_grid(vec![Vec3::ZERO; 2], vec![0.0, 180.0]);
1806        assert!(
1807            root_trajectory_metrics(&half_turn_step, 0)
1808                .unwrap()
1809                .yaw
1810                .is_none()
1811        );
1812    }
1813
1814    #[test]
1815    fn root_yaw_uses_animated_ancestor_model_rotation_for_fixed_child_local_pose() {
1816        let skeleton = Skeleton {
1817            bones: vec![
1818                Bone {
1819                    name: "ancestor".into(),
1820                    parent: None,
1821                    rest: Transform::IDENTITY,
1822                    inverse_bind: None,
1823                },
1824                Bone {
1825                    name: "root".into(),
1826                    parent: Some(0),
1827                    rest: Transform::IDENTITY,
1828                    inverse_bind: None,
1829                },
1830            ],
1831        };
1832        let clip = Clip {
1833            name: "inherited_yaw".into(),
1834            duration_s: 1.0,
1835            tracks: vec![
1836                Track {
1837                    bone: 0,
1838                    property: Property::Rotation,
1839                    interpolation: Interpolation::Linear,
1840                    times: vec![0.0, 0.5, 1.0],
1841                    values: TrackValues::Quats(
1842                        [0.0_f32, 45.0, 90.0]
1843                            .map(|degrees| Quat::from_rotation_y(degrees.to_radians()))
1844                            .to_vec(),
1845                    ),
1846                },
1847                Track {
1848                    bone: 1,
1849                    property: Property::Rotation,
1850                    interpolation: Interpolation::Linear,
1851                    times: vec![0.0, 0.5, 1.0],
1852                    values: TrackValues::Quats(vec![Quat::IDENTITY; 3]),
1853                },
1854            ],
1855        };
1856        let grid = sample_clip(&skeleton, &clip, 3);
1857        for frame in 0..grid.frame_count() {
1858            assert_eq!(grid.local(frame, 1).rotation, Quat::IDENTITY);
1859        }
1860
1861        let trajectory = root_trajectory_metrics(&grid, 1).expect("child trajectory");
1862        assert_stationary_translation(&trajectory);
1863        let yaw = trajectory.yaw.expect("ancestor supplies model-space yaw");
1864        assert_eq!(yaw.heading_axis, RootYawHeadingAxis::PositiveZ);
1865        assert!((yaw.net_yaw_deg - 90.0).abs() < 1.0e-4);
1866        assert!((yaw.unwrapped_yaw_deg - 90.0).abs() < 1.0e-4);
1867        assert!((yaw.yaw_travel_deg - 90.0).abs() < 1.0e-4);
1868    }
1869}