完善Phase-A会话级联合优化并修正雷达相位中心高度先验

Co-authored-by: Cursor <cursoragent@cursor.com>

README.md

@@ -44,12 +44,15 @@ p_IMU = T_IMU_lidar · p_lidar
 | `summary.json`     | 状态、残差、可观性              |
 
 
-新车原始数据导出（H32 dlog + N300 rscap）：
+新车原始数据导出（H32 dlog/zip + HI13 rscap）：
 
 ```powershell
 python tools\export_rscap_to_v1.py `
-  --imu-rscap path\to\n300.rscap `
-  --lidar-dlog path\to\session_or_dlog `
+  --imu-rscap path\to\hi13r4-imu.rscap `
+  --imu-kind hi13 `
+  --lidar-dlog path\to\session_or_recovered.zip `
+  --host-start 2026-08-08T17:40:05 `
+  --host-end 2026-08-08T17:45:15 `
   --out path\to\session_v1 `
   --require-difop
 ```
@@ -70,16 +73,25 @@ powershell -File tools\reproduce_synthetic.ps1
 证明：链路可跑通，能收回已知 yaw / δt。  
 不证明：实车安装精度、平移可交付。
 
-产物在 `examples/synthetic_session/out/`。叠点查看：
+产物在 `examples/synthetic_session/out/`（含 `summary.json` 与 `motion_pairs.json`）。叠点查看：
 
 ```powershell
+# 优先读取 summary 同目录的 motion_pairs.json，按需加载点云（无需重算配准）
 python tools\visualize_pair_3d.py `
   --lidar examples\synthetic_session\lidar `
-  --imu examples\synthetic_session\imu.csv `
   --summary examples\synthetic_session\out\summary.json `
   --pair-index 0
 ```
 
+旧标定目录若缺少缓存，可只补导出运动对（不重求解外参）：
+
+```powershell
+python tools\export_motion_pairs_for_viz.py `
+  --lidar path\to\lidar `
+  --imu path\to\imu.csv `
+  --summary path\to\out\summary.json
+```
+
 键 `1`–`4` 切换叠点模式；`N`/`P` 切换运动对。
 
 ---

config/vehicle_hi13_h32_20260808.yaml

@@ -0,0 +1,86 @@
+schema_version: 1
+
+vehicle:
+  vehicle_id: "outdoor_usable_20260808"
+  body_frame:
+    name: "base_link"
+    # 车体约定：后轮轴中心在地面投影为原点附近参考；X 前 / Y 左 / Z 上
+    # translation_m 的 Z 使用「离地高度」；后轮轴中心离地 294 mm
+    axes: "X forward, Y left, Z up"
+    unit: m
+    reference_point: "rear_axle_center_xy__z_above_ground"
+    rear_axle_height_above_ground_m: 0.294
+
+installation:
+  installation_id: "20260808_priority_windows"
+  installed_at: "2026-08-08"
+  notes: >
+    HI13R4 + H32 DLogCapture. Body +X forward: LiDAR and IMU at positive X.
+    CAD sheet may draw +X rearward; numbers below are body-frame.
+    Z is height above ground = CAD height at axle + 0.294 m (axle AGL).
+    LiDAR CAD dZ is 1637.499879 mm relative to the axle reference. Phase-center
+    AGL adds rear-axle height 294 mm and the 63.5 mm phase-center offset.
+    IMU axes: HI13R4 manual §2.4 RFU (X right, Y forward, Z up).
+    LiDAR Cartesian in NPZ assumed body-aligned (X forward).
+
+sensors:
+  imu:
+    model: "HI13R4"
+    raw_frame:
+      # HI13R4 用户手册 2.4：右-前-上 (RFU)
+      axes: "X right, Y forward, Z up (RFU)"
+      driver_axis_remapped: false
+    mount_in_body:
+      # X/Y：后轮轴中心 → IMU；Z：离地 = CAD 0.8925 + 0.294
+      translation_m: [2.574126255, 0.0365, 1.1865]
+      # body <- imu : p_body = R_body_imu * p_imu
+      # R_body_imu = [[0,1,0],[-1,0,0],[0,0,1]]  (fwd=imu_y, left=-imu_x, up=imu_z)
+      rotation_matrix_body_imu: [[0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]
+      rotation_quaternion_xyzw: null
+      source: "CAD X/Y in body (+X forward); Z = CAD axle-height + 294mm AGL + HI13R4 RFU"
+
+  lidar:
+    model: "RSLidarH32"
+    points_field: points
+    raw_frame:
+      axes: "X forward, Y left, Z up (Cartesian metres in NPZ points)"
+      driver_axis_remapped: false
+    mount_in_body:
+      # X/Y：后轮轴中心 → 雷达
+      # Z离地 = CAD dZ 1.637499879 + 后轮轴离地 0.294 + 相位中心偏移 0.0635
+      translation_m: [2.522276859, 0.000020526, 1.994999879]
+      rotation_matrix_body_lidar: [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
+      rotation_quaternion_xyzw: null
+      source: "CAD X/Y in body (+X forward); Z AGL = CAD dZ 1.637499879 + axle AGL 0.294 + phase-center offset 0.0635; attitude = body"
+
+  rtk:
+    frame_definition: ""
+    reference_point: ""
+    existing_T_RTK_LIDAR_file: ""
+
+time:
+  imu_timestamp_source: "hi13_device_timestamp_ms_seconds"
+  lidar_timestamp_source: "h32_msop_device_timestamp_seconds"
+  lidar_frame_time_definition: "t_start/t_end in frames_index.csv; pipeline uses midpoint"
+  host_bridge: "MSOP HostReceiveUtcTicks + IMU receive_utc_ticks"
+
+# Derived prior for p_IMU = R_IMU_lidar * p_lidar + t_IMU_lidar
+# t_body = t_lidar_body - t_imu_body
+# t_IMU_lidar = R_IMU_body * t_body,  R_IMU_lidar = R_IMU_body * R_body_lidar
+derived_T_IMU_lidar_prior:
+  R_IMU_lidar: [[0.0, -1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]
+  t_IMU_lidar_m: [0.036479474, -0.051849396, 0.808499879]
+  t_lidar_from_imu_in_body_m: [-0.051849396, -0.036479474, 0.808499879]
+  notes: >
+    Rotation prior ~90 deg yaw (body/lidar X-fwd vs IMU Y-fwd).
+    Relative Z = 1.994999879 - 1.1865 = 0.808499879 m.
+
+initialization:
+  translation_prior:
+    enabled: true
+    sigma_m: [0.05, 0.05, 0.05]
+    t_IMU_lidar_m: [0.036479474, -0.051849396, 0.808499879]
+  rotation_prior:
+    enabled: true
+    sigma_deg: 15.0
+    R_IMU_lidar: [[0.0, -1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]

docs/20260808_LiDAR-IMU标定现状与问题.md

@@ -0,0 +1,209 @@
+# 20260808 HI13 + H32：LiDAR–IMU 标定现状与问题
+
+> 数据：`D:\data\calibration_usable_20260808`  
+> 可用会话：`sessions_v1_host_aligned`（三优先窗）  
+> 当前结果目录：各窗 `out_fixed_dt0/`  
+> 清单：`sessions_v1_host_aligned/calibration_manifest_fixed_dt0.json`  
+> 车辆配置：`config/vehicle_hi13_h32_20260808.yaml`  
+> 约定外参：`p_IMU = T_IMU_lidar · p_lidar`
+
+---
+
+## 1. 一句话结论
+
+**旋转 + 主机桥接时间对齐可以冻结；平移（full_se3）尚不可正式交付。**  
+三窗 `rotation_only`（δt=0）结果跨窗一致，**不必因平移先验 Z 修正而重跑旋转**。
+
+---
+
+## 2. 当前可用结果（`out_fixed_dt0`）
+
+约定：`p_IMU = T_IMU_lidar · p_lidar`；本轮交付 **仅旋转**，`t = [0,0,0]`，`time_offset_s = 0`。
+
+| 窗 | 状态 | δt | roll/pitch/yaw (°) | 手眼 RMS (°) | 手眼对数 | vs CAD prior |
+|----|------|----|---------------------|--------------|----------|--------------|
+| `priority_174005_174515` | `rotation_only_accepted` | 0 | −0.398 / +0.108 / **89.998** | 0.625 | 1374 | 0.413° |
+| `priority_174905_175450` | 同上 | 0 | −0.316 / −0.352 / **90.002** | 0.293 | 1182 | 0.473° |
+| `priority_175910_180530` | 同上 | 0 | −0.373 / −0.036 / **90.005** | 0.786 | 789 | 0.375° |
+
+- 跨窗旋转互差约 **0.15°–0.47°**（相对三窗均值 ≤0.26°）。
+- CAD/安装平移先验只用于后续 SE3 / 校验，不写入本轮交付 `T`。
+- 原始摘要：各窗 `out_fixed_dt0/summary.json`；总表 `calibration_manifest_fixed_dt0.json`。
+
+### 2.1 窗1 `priority_174005_174515` — `R_IMU_lidar`
+
+- 路径：`...\priority_174005_174515\out_fixed_dt0\summary.json`
+- rpy_deg_xyz：`[-0.39806616272552936, 0.10842366761721789, 89.99788311264182]`
+- quaternion_xyzw：`[-0.0031254048203223084, -0.0017872288323561246, 0.7070914597978358, 0.707112936622415]`
+
+```text
+R =
+[[ 3.6946588133e-05, -0.9999758655693408, -0.0069474393698490 ],
+ [ 0.9999982088237712,  2.3798651350e-05,  0.0018925598731371 ],
+ [-0.0018923488575947, -0.0069474968493909,  0.9999740753156198 ]]
+t = [0, 0, 0]
+```
+
+### 2.2 窗2 `priority_174905_175450` — `R_IMU_lidar`
+
+- 路径：`...\priority_174905_175450\out_fixed_dt0\summary.json`
+- rpy_deg_xyz：`[-0.31554362742542746, -0.35231680831805484, 90.00193338170823]`
+- quaternion_xyzw：`[0.00022698471903919405, -0.004121119694188399, 0.7071067019847445, 0.7070948146172911]`
+
+```text
+R =
+[[-3.3743238552e-05, -0.9999848355714863, -0.0055070399001942 ],
+ [ 0.9999810938467022,  1.2097239027e-07, -0.0061491421465438 ],
+ [ 0.0061490495645171, -0.0055071432752239,  0.9999659297008071 ]]
+t = [0, 0, 0]
+```
+
+### 2.3 窗3 `priority_175910_180530` — `R_IMU_lidar`
+
+- 路径：`...\priority_175910_180530\out_fixed_dt0\summary.json`
+- rpy_deg_xyz：`[-0.37335575043737196, -0.03585856981709423, 90.00538974486676]`
+- quaternion_xyzw：`[-0.002082462199099642, -0.002525219418619018, 0.7071355299053291, 0.7070704554452735]`
+
+```text
+R =
+[[-9.4068775206e-05, -0.9999787650354249, -0.0065161821101807 ],
+ [ 0.9999997997313592, -8.9988606603e-05, -0.0006264497522950 ],
+ [ 0.0006258500675082, -0.0065162397345547,  0.9999785732361544 ]]
+t = [0, 0, 0]
+```
+
+### 相对历史失败轮次
+
+| 轮次 | 问题 | 结果 |
+|------|------|------|
+| `sessions_v1_aligned` | 首帧强行对齐设备钟 | 三窗手眼失败，RMS ~9°–12° |
+| 自由估 δt + signed refine | 窗3 δt 漂到 −0.48 s；窗2 yaw≈19° | 跨窗 yaw 矛盾（81°/19°/93°） |
+| **本轮 fixed δt=0** | 主机桥接后冻结时间 | 三窗 yaw≈90°，可互证 |
+
+---
+
+## 3. 已澄清并写入配置的坐标系 / 先验
+
+### 3.1 车体与传感器
+
+- 车体：X 前 / Y 左 / Z 上；雷达与 IMU 安装在 **X 正方向**（后轮轴前方）。
+- CAD 图纸可能画成 +X 朝后，那只是读图坐标系，**不是**车体真实轴。
+- IMU：HI13 RFU（X 右 / Y 前 / Z 上），原始数据不做轴向重映射。
+- 雷达 NPZ：假定与车体一致（X 前 / Y 左 / Z 上）。
+
+### 3.2 安装量（`translation_m`）
+
+| 传感器 | X / Y（后轮轴中心） | Z（离地） |
+|--------|---------------------|-----------|
+| IMU | 2.574 / 0.0365 m | 0.8925 + 0.294 = **1.1865 m** |
+| 雷达 | 2.522 / 0.00002 m | 相位中心离地 **1.994999879 m** |
+
+- 后轮轴中心离地：**294 mm**（Z 用离地高时加在 CAD 轴心高上）。
+- 雷达 CAD `dZ=1.637499879 m`；相位中心离地还需加后轮轴离地 `0.294 m` 和相位中心偏移 `0.0635 m`，最终为 `1.994999879 m`。
+
+### 3.3 导出外参先验
+
+- `R_IMU_lidar` ≈ yaw 90°：`[[0,-1,0],[1,0,0],[0,0,1]]`（软约束 σ=15°）。
+- `t_IMU_lidar` ≈ **`[0.0365, -0.0518, 0.8085]` m**（相对 Z = 1.994999879 − 1.1865 = 0.808499879 m）。
+- **旋转先验不因 Z 修正改变**；平移先验 Z 更新为 0.808499879 m。
+
+---
+
+## 4. 现存问题清单
+
+### P1. IMU 预积分平移 `Δp` 不可用（阻塞正式平移）
+
+- 现象：可视化模式 4 若用完整 `X⁻¹ A X`，橙/蓝点云常呈**上下错层**（Z 差米级～几十米）。
+- 根因：加速度预积分缺少可靠重力/零偏处理，`t_A` 尤其 Z 发散；**不是旋转外参错了**。
+- 旁证：相对 GICP 的旋转残差中位约 0.16°；`|t_A|` 中位却常 >1 m。
+- 影响：`full_se3` / 依赖 IMU 位移的平移估计不可信。
+- 缓解（已做）：`visualize_pair_3d.py` 对 `rotation_only` 默认模式 4 = **R 共轭 + GICP 的 t_B**（`--mode4-translation gicp|imu|auto`）。
+
+### P2. 平面运动导致竖直平移弱可观
+
+- 三优先窗以水平转弯为主，缺少缓坡/俯仰激励。
+- 流水线门控已给出 `translation_accepted=false`。
+- 即使打开平移先验（σ≈5 cm），弱激励下结果易变成**先验回显**，不宜当标定成功。
+
+### P3. 时间偏移若再自由估计会被带偏（已规避，需保持）
+
+- 主机 UTC 桥接（MSOP/IMU `HostReceiveUtc`）后，两路已在同一时间轴，残差通常几十毫秒量级。
+- 若再做有符号 δt 精修，会与错误/未收敛的 R 耦合，窗3 曾从约 −0.12 s 走到 **−0.48 s**。
+- **现行做法**：桥接会话使用 `--fixed-time-offset-s 0 --no-signed-time-refine`。
+
+### P4. 单窗低残差 ≠ 外参正确（历史教训）
+
+- 自由 δt 轮次中，窗2 手眼 RMS 最低（~0.3°）但 yaw≈19°，与 CAD/其他窗差 60°+。
+- 平面运动下 yaw 外参可出现多个能拟合 `R_A R_X ≈ R_X R_B` 的解。
+- **必须**做跨窗一致性 + 可视化叠点，不能只看单窗 RMS。
+
+### P5. 旋转软先验尚未做无先验对照
+
+- 当前 σ=15°；笔记显示 Tsai 初值本身已接近（约 0.3°–1.1° RMS），不像纯先验硬拽。
+- 仍缺一次：关闭先验或放大 `sigma_deg` 的对照，以排除「只是被拉到 90°」的疑虑。
+
+### P6. 文档与操作约定未完全同步（工程）
+
+- README 需明确写清：host-bridge 后固定 δt=0、禁用 signed refine、rotation_only 可视化用法。
+- 交付物目前缺一版「冻结的联合/中位 R + 使用说明」JSON/报告（旋转可交，平移明确不交）。
+
+---
+
+## 5. 不该做 / 可以做
+
+| 动作 | 建议 |
+|------|------|
+| 因 Z 先验修正重跑三窗 rotation_only | **不必**（R 未依赖新 t） |
+| 正式交付 6-DOF / 信赖当前 `Δp` 估 t | **不要** |
+| 试验性 `full_se3`（固定 R、δt=0、新 t 先验） | 可做，结果标「实验」 |
+| 可视化验收模式 3 vs 4（gicp 平移） | **建议做** |
+| 无先验 / 大 σ 旋转对照 | **建议做** |
+| 冻结交付 `R` + `δt=0` 说明 | **建议做** |
+| 补采缓坡或加强垂直尺寸约束后再估 t | 正式平移前需要 |
+
+---
+
+## 6. 建议下一步顺序
+
+1. **验收旋转**：三窗抽转弯运动对，模式 3/4 叠点；可选无先验对照。  
+2. **定稿旋转**：三窗中位或联合手眼 → 交付 `R_IMU_lidar` +「δt=0（主机桥接）」说明；**明确不交 t**。  
+3. **工程收尾**：README 主机桥接配方；需要时再整理联合标定脚本入口。  
+4. **平移（靠后）**：改善 IMU 位移模型或改用更可靠的位移观测 + 竖直激励后，再用新 `t` 先验跑 SE3。
+
+---
+
+## 7. 常用路径与命令
+
+```text
+数据根：
+  D:\data\calibration_usable_20260808\sessions_v1_host_aligned\
+
+结果：
+  ...\priority_XXXX\out_fixed_dt0\summary.json
+  ...\priority_XXXX\out_fixed_dt0\motion_pairs.json
+  ...\calibration_manifest_fixed_dt0.json
+```
+
+```powershell
+# 可视化（rotation_only 默认模式4用 GICP 平移）
+python tools\visualize_pair_3d.py `
+  --lidar D:\data\calibration_usable_20260808\sessions_v1_host_aligned\priority_174005_174515\lidar `
+  --summary D:\data\calibration_usable_20260808\sessions_v1_host_aligned\priority_174005_174515\out_fixed_dt0\summary.json `
+  --pair-index 0
+
+# 若要看「坏 Δp」导致的错层效果：
+#   --mode4-translation imu
+```
+
+---
+
+## 8. 问题优先级（跟踪用）
+
+| ID | 严重度 | 状态 | 标题 |
+|----|--------|------|------|
+| P1 | 高 | 未解决 | IMU `Δp` 不可用，阻塞正式平移 |
+| P2 | 高 | 未解决 | 平面运动，竖直 t 弱可观 |
+| P3 | 高 | 已规避 | 自由 δt / signed refine 带偏（需保持冻结） |
+| P4 | 中 | 已吸收教训 | 单窗低残差不可单独验收 |
+| P5 | 中 | 待做 | 无旋转先验对照 |
+| P6 | 低 | 待做 | README/交付物同步 |

docs/V1_数据格式.md

@@ -4,20 +4,24 @@
 
 ## 从原始数据导出
 
-**推荐（新 H32 插件 `RSLidarH32_3D_DLogCaptureNet48`）：** N300 `.rscap` + 雷达 Medulla dlog（含 raw MSOP / DIFOP）。
+**推荐（新 H32 + HI13）：** HI13 `.rscap` + 雷达 Medulla dlog / recovered zip（raw MSOP + DIFOP）。
 
 ```powershell
 python tools\export_rscap_to_v1.py `
-  --imu-rscap path\to\n300.rscap `
-  --lidar-dlog path\to\session_or_dlog `
+  --imu-rscap path\to\hi13r4-imu.rscap `
+  --imu-kind hi13 `
+  --lidar-dlog path\to\session_or_dlog_or_recovered.zip `
+  --host-start 2026-08-08T17:40:05 `
+  --host-end 2026-08-08T17:45:15 `
   --out path\to\session_v1 `
-  --frame-stride 1 `
+  --frame-stride 5 `
   --require-difop
 ```
 
-`--lidar-dlog` 指向含 `dobject/` + `dobject_recording/` 的目录（或其上级含 `dlog/` 子目录亦可）。  
-默认 DObject：`frontlidar-msop-raw`、`frontlidar-difop-raw`（可用 `--msop-object` / `--difop-object` 覆盖）。  
-有 DIFOP 时用设备通道角做 XYZ；`--require-difop` 在缺少有效 DIFOP 时直接失败。
+`--lidar-dlog` 可为：标准 `dobject/`+`dobject_recording/` 目录，或 recovered zip（`indices.log` + `data.bin`）。  
+`--imu-kind`：`hi13` / `n300` / `auto`（默认按文件名推断）。  
+`--host-start/end`：按本地墙钟切窗（仅裁剪；标定主轴仍是设备时间）。  
+默认 DObject：`frontlidar-msop-raw`、`frontlidar-difop-raw`。
 
 **兼容旧 MSOP-only `.rscap`：**
 
@@ -30,7 +34,7 @@ python tools\export_rscap_to_v1.py `
 ```
 
 产出：`imu.csv`、`lidar/`（含 `frames_index.csv`）、`export_summary.json`。  
-时间轴为**设备时间**：N300 `device_timestamp_us`→秒；H32 MSOP 设备时间戳→秒。主机接收时间不写入标定主轴。
+标定主轴仍是**设备时间**；同时写出**主机 UTC 接收时间**，用于把雷达帧桥接到 IMU 设备钟（禁止把两边设备时间第一帧强行重合）。
 
 ## IMU
 
@@ -39,16 +43,18 @@ python tools\export_rscap_to_v1.py `
 ### CSV
 
 ```text
-t,gx,gy,gz,ax,ay,az
-0.000000000,0.01,-0.02,0.00,0.05,-0.03,9.81
+t,gx,gy,gz,ax,ay,az,t_host_utc_s,receive_utc_ticks
+0.000000000,0.01,-0.02,0.00,0.05,-0.03,9.81,1754646005.123,6389...
 ...
 ```
 
 | 列 | 含义 | 单位 |
 |---|---|---|
-| t | IMU 时钟时间 | s |
+| t | IMU 设备时钟时间 | s |
 | gx,gy,gz | 角速度 | rad/s |
 | ax,ay,az | 比力/加速度 | m/s² |
+| t_host_utc_s | 主机 UTC 接收时间（Unix） | s |
+| receive_utc_ticks | 同上，.NET UTC ticks | — |
 
 ### NPZ
 
@@ -72,12 +78,16 @@ lidar_session/
 ### frames_index.csv
 
 ```text
-frame_id,filename,t_start,t_end
-0,frames/frame_00000.npz,10.000,10.100
-1,frames/frame_00001.npz,10.100,10.200
+frame_id,filename,t_start,t_end,host_receive_utc_ticks,t_host_utc_s,host_receive_utc_end_ticks,t_host_utc_end_s
+0,frames/frame_00000.npz,10.000,10.100,6389...,1754646005.12,6389...,1754646005.22
 ```
 
-也兼容旧列名 `file`（NumPy 读取时可能变成 `file_`）。
+| 列 | 含义 |
+|---|---|
+| t_start / t_end | H32 MSOP **设备时间**（秒） |
+| t_host_utc_s / t_host_utc_end_s | 帧首/末包 **HostReceiveUtcTicks** → Unix 秒 |
+
+也兼容旧列名 `file`。对齐脚本用主机 UTC 把 `t_*` 重写到 IMU 设备钟后再跑标定。
 
 ### 每帧 NPZ
 

docs/问题清单.md

@@ -0,0 +1,209 @@
+# 20260808 HI13 + H32：LiDAR–IMU 标定现状与问题
+
+> 数据：`D:\data\calibration_usable_20260808`  
+> 可用会话：`sessions_v1_host_aligned`（三优先窗）  
+> 当前结果目录：各窗 `out_fixed_dt0/`  
+> 清单：`sessions_v1_host_aligned/calibration_manifest_fixed_dt0.json`  
+> 车辆配置：`config/vehicle_hi13_h32_20260808.yaml`  
+> 约定外参：`p_IMU = T_IMU_lidar · p_lidar`
+
+---
+
+## 1. 一句话结论
+
+**旋转 + 主机桥接时间对齐可以冻结；平移（full_se3）尚不可正式交付。**  
+三窗 `rotation_only`（δt=0）结果跨窗一致，**不必因平移先验 Z 修正而重跑旋转**。
+
+---
+
+## 2. 当前可用结果（`out_fixed_dt0`）
+
+约定：`p_IMU = T_IMU_lidar · p_lidar`；本轮交付 **仅旋转**，`t = [0,0,0]`，`time_offset_s = 0`。
+
+| 窗 | 状态 | δt | roll/pitch/yaw (°) | 手眼 RMS (°) | 手眼对数 | vs CAD prior |
+|----|------|----|---------------------|--------------|----------|--------------|
+| `priority_174005_174515` | `rotation_only_accepted` | 0 | −0.398 / +0.108 / **89.998** | 0.625 | 1374 | 0.413° |
+| `priority_174905_175450` | 同上 | 0 | −0.316 / −0.352 / **90.002** | 0.293 | 1182 | 0.473° |
+| `priority_175910_180530` | 同上 | 0 | −0.373 / −0.036 / **90.005** | 0.786 | 789 | 0.375° |
+
+- 跨窗旋转互差约 **0.15°–0.47°**（相对三窗均值 ≤0.26°）。
+- CAD/安装平移先验只用于后续 SE3 / 校验，不写入本轮交付 `T`。
+- 原始摘要：各窗 `out_fixed_dt0/summary.json`；总表 `calibration_manifest_fixed_dt0.json`。
+
+### 2.1 窗1 `priority_174005_174515` — `R_IMU_lidar`
+
+- 路径：`...\priority_174005_174515\out_fixed_dt0\summary.json`
+- rpy_deg_xyz：`[-0.39806616272552936, 0.10842366761721789, 89.99788311264182]`
+- quaternion_xyzw：`[-0.0031254048203223084, -0.0017872288323561246, 0.7070914597978358, 0.707112936622415]`
+
+```text
+R =
+[[ 3.6946588133e-05, -0.9999758655693408, -0.0069474393698490 ],
+ [ 0.9999982088237712,  2.3798651350e-05,  0.0018925598731371 ],
+ [-0.0018923488575947, -0.0069474968493909,  0.9999740753156198 ]]
+t = [0, 0, 0]
+```
+
+### 2.2 窗2 `priority_174905_175450` — `R_IMU_lidar`
+
+- 路径：`...\priority_174905_175450\out_fixed_dt0\summary.json`
+- rpy_deg_xyz：`[-0.31554362742542746, -0.35231680831805484, 90.00193338170823]`
+- quaternion_xyzw：`[0.00022698471903919405, -0.004121119694188399, 0.7071067019847445, 0.7070948146172911]`
+
+```text
+R =
+[[-3.3743238552e-05, -0.9999848355714863, -0.0055070399001942 ],
+ [ 0.9999810938467022,  1.2097239027e-07, -0.0061491421465438 ],
+ [ 0.0061490495645171, -0.0055071432752239,  0.9999659297008071 ]]
+t = [0, 0, 0]
+```
+
+### 2.3 窗3 `priority_175910_180530` — `R_IMU_lidar`
+
+- 路径：`...\priority_175910_180530\out_fixed_dt0\summary.json`
+- rpy_deg_xyz：`[-0.37335575043737196, -0.03585856981709423, 90.00538974486676]`
+- quaternion_xyzw：`[-0.002082462199099642, -0.002525219418619018, 0.7071355299053291, 0.7070704554452735]`
+
+```text
+R =
+[[-9.4068775206e-05, -0.9999787650354249, -0.0065161821101807 ],
+ [ 0.9999997997313592, -8.9988606603e-05, -0.0006264497522950 ],
+ [ 0.0006258500675082, -0.0065162397345547,  0.9999785732361544 ]]
+t = [0, 0, 0]
+```
+
+### 相对历史失败轮次
+
+| 轮次 | 问题 | 结果 |
+|------|------|------|
+| `sessions_v1_aligned` | 首帧强行对齐设备钟 | 三窗手眼失败，RMS ~9°–12° |
+| 自由估 δt + signed refine | 窗3 δt 漂到 −0.48 s；窗2 yaw≈19° | 跨窗 yaw 矛盾（81°/19°/93°） |
+| **本轮 fixed δt=0** | 主机桥接后冻结时间 | 三窗 yaw≈90°，可互证 |
+
+---
+
+## 3. 已澄清并写入配置的坐标系 / 先验
+
+### 3.1 车体与传感器
+
+- 车体：X 前 / Y 左 / Z 上；雷达与 IMU 安装在 **X 正方向**（后轮轴前方）。
+- CAD 图纸可能画成 +X 朝后，那只是读图坐标系，**不是**车体真实轴。
+- IMU：HI13 RFU（X 右 / Y 前 / Z 上），原始数据不做轴向重映射。
+- 雷达 NPZ：假定与车体一致（X 前 / Y 左 / Z 上）。
+
+### 3.2 安装量（`translation_m`）
+
+| 传感器 | X / Y（后轮轴中心） | Z（离地） |
+|--------|---------------------|-----------|
+| IMU | 2.574 / 0.0365 m | 0.8925 + 0.294 = **1.1865 m** |
+| 雷达 | 2.522 / 0.00002 m | 相位中心离地 **1.994999879 m** |
+
+- 后轮轴中心离地：**294 mm**（Z 用离地高时加在 CAD 轴心高上）。
+- 雷达 CAD `dZ=1.637499879 m`；相位中心离地还需加后轮轴离地 `0.294 m` 和相位中心偏移 `0.0635 m`，最终为 `1.994999879 m`。
+
+### 3.3 导出外参先验
+
+- `R_IMU_lidar` ≈ yaw 90°：`[[0,-1,0],[1,0,0],[0,0,1]]`（软约束 σ=15°）。
+- `t_IMU_lidar` ≈ **`[0.0365, -0.0518, 0.8085]` m**（相对 Z = 1.994999879 − 1.1865 = 0.808499879 m）。
+- **旋转先验不因 Z 修正改变**；平移先验 Z 更新为 0.808499879 m。
+
+---
+
+## 4. 现存问题清单
+
+### P1. IMU 预积分平移 `Δp` 不可用（阻塞正式平移）
+
+- 现象：可视化模式 4 若用完整 `X⁻¹ A X`，橙/蓝点云常呈**上下错层**（Z 差米级～几十米）。
+- 根因：加速度预积分缺少可靠重力/零偏处理，`t_A` 尤其 Z 发散；**不是旋转外参错了**。
+- 旁证：相对 GICP 的旋转残差中位约 0.16°；`|t_A|` 中位却常 >1 m。
+- 影响：`full_se3` / 依赖 IMU 位移的平移估计不可信。
+- 缓解（已做）：`visualize_pair_3d.py` 对 `rotation_only` 默认模式 4 = **R 共轭 + GICP 的 t_B**（`--mode4-translation gicp|imu|auto`）。
+
+### P2. 平面运动导致竖直平移弱可观
+
+- 三优先窗以水平转弯为主，缺少缓坡/俯仰激励。
+- 流水线门控已给出 `translation_accepted=false`。
+- 即使打开平移先验（σ≈5 cm），弱激励下结果易变成**先验回显**，不宜当标定成功。
+
+### P3. 时间偏移若再自由估计会被带偏（已规避，需保持）
+
+- 主机 UTC 桥接（MSOP/IMU `HostReceiveUtc`）后，两路已在同一时间轴，残差通常几十毫秒量级。
+- 若再做有符号 δt 精修，会与错误/未收敛的 R 耦合，窗3 曾从约 −0.12 s 走到 **−0.48 s**。
+- **现行做法**：桥接会话使用 `--fixed-time-offset-s 0 --no-signed-time-refine`。
+
+### P4. 单窗低残差 ≠ 外参正确（历史教训）
+
+- 自由 δt 轮次中，窗2 手眼 RMS 最低（~0.3°）但 yaw≈19°，与 CAD/其他窗差 60°+。
+- 平面运动下 yaw 外参可出现多个能拟合 `R_A R_X ≈ R_X R_B` 的解。
+- **必须**做跨窗一致性 + 可视化叠点，不能只看单窗 RMS。
+
+### P5. 旋转软先验尚未做无先验对照
+
+- 当前 σ=15°；笔记显示 Tsai 初值本身已接近（约 0.3°–1.1° RMS），不像纯先验硬拽。
+- 仍缺一次：关闭先验或放大 `sigma_deg` 的对照，以排除「只是被拉到 90°」的疑虑。
+
+### P6. 文档与操作约定未完全同步（工程）
+
+- README 需明确写清：host-bridge 后固定 δt=0、禁用 signed refine、rotation_only 可视化用法。
+- 交付物目前缺一版「冻结的联合/中位 R + 使用说明」JSON/报告（旋转可交，平移明确不交）。
+
+---
+
+## 5. 不该做 / 可以做
+
+| 动作 | 建议 |
+|------|------|
+| 因 Z 先验修正重跑三窗 rotation_only | **不必**（R 未依赖新 t） |
+| 正式交付 6-DOF / 信赖当前 `Δp` 估 t | **不要** |
+| 试验性 `full_se3`（固定 R、δt=0、新 t 先验） | 可做，结果标「实验」 |
+| 可视化验收模式 3 vs 4（gicp 平移） | **建议做** |
+| 无先验 / 大 σ 旋转对照 | **建议做** |
+| 冻结交付 `R` + `δt=0` 说明 | **建议做** |
+| 补采缓坡或加强垂直尺寸约束后再估 t | 正式平移前需要 |
+
+---
+
+## 6. 建议下一步顺序
+
+1. **验收旋转**：三窗抽转弯运动对，模式 3/4 叠点；可选无先验对照。  
+2. **定稿旋转**：三窗中位或联合手眼 → 交付 `R_IMU_lidar` +「δt=0（主机桥接）」说明；**明确不交 t**。  
+3. **工程收尾**：README 主机桥接配方；需要时再整理联合标定脚本入口。  
+4. **平移（靠后）**：改善 IMU 位移模型或改用更可靠的位移观测 + 竖直激励后，再用新 `t` 先验跑 SE3。
+
+---
+
+## 7. 常用路径与命令
+
+```text
+数据根：
+  D:\data\calibration_usable_20260808\sessions_v1_host_aligned\
+
+结果：
+  ...\priority_XXXX\out_fixed_dt0\summary.json
+  ...\priority_XXXX\out_fixed_dt0\motion_pairs.json
+  ...\calibration_manifest_fixed_dt0.json
+```
+
+```powershell
+# 可视化（rotation_only 默认模式4用 GICP 平移）
+python tools\visualize_pair_3d.py `
+  --lidar D:\data\calibration_usable_20260808\sessions_v1_host_aligned\priority_174005_174515\lidar `
+  --summary D:\data\calibration_usable_20260808\sessions_v1_host_aligned\priority_174005_174515\out_fixed_dt0\summary.json `
+  --pair-index 0
+
+# 若要看「坏 Δp」导致的错层效果：
+#   --mode4-translation imu
+```
+
+---
+
+## 8. 问题优先级（跟踪用）
+
+| ID | 严重度 | 状态 | 标题 |
+|----|--------|------|------|
+| P1 | 高 | 未解决 | IMU `Δp` 不可用，阻塞正式平移 |
+| P2 | 高 | 未解决 | 平面运动，竖直 t 弱可观 |
+| P3 | 高 | 已规避 | 自由 δt / signed refine 带偏（需保持冻结） |
+| P4 | 中 | 已吸收教训 | 单窗低残差不可单独验收 |
+| P5 | 中 | 待做 | 无旋转先验对照 |
+| P6 | 低 | 待做 | README/交付物同步 |

imu_lidar/CHANGELOG.md

@@ -5,6 +5,45 @@
 
 ---
 
+## 2026-08-11 10:55 (UTC+8)
+
+### 运动对缓存：标定落盘，可视化直读
+
+- **原本**：`visualize_pair_3d` 每次启动都重新关键帧+配准+预积分，等同半次标定。
+- **改成**：
+  - 标定成功后写出 `motion_pairs.json`（`motion_pairs_io.py` / `finalize`）。
+  - 可视化优先读缓存并对点云懒加载；`--rebuild-pairs` 可回退旧路径。
+  - 旧结果可用 `tools/export_motion_pairs_for_viz.py` 只补导出运动对，无需重求解外参。
+
+---
+
+## 2026-08-11 08:55 (UTC+8)
+
+### 主机桥接后冻结 δt + 旋转先验软约束
+
+- **原本**：手眼后 signed δt 精修可在弱 MSE 下降下连走数步（最远约 0.5 s）；旋转手眼无 CAD 先验，平面运动下 yaw 易掉进低残差错解。
+- **改成**：
+  - CLI：`--fixed-time-offset-s`、`--no-signed-time-refine`、`--max-signed-refine-shift-s`。
+  - signed refine：默认 `|Δδt|≤0.05 s`，且要求 MSE 至少降约 2%。
+  - `rotation_handeye` 读取配置 `rotation_prior` 作初值/软约束。
+  - 主机 UTC 桥接会话建议：`--fixed-time-offset-s 0 --no-signed-time-refine`。
+
+---
+
+## 2026-08-09 14:30 (UTC+8)
+
+### 导出：HI13 IMU + recovered dlog zip + 墙钟切窗
+
+- **原本**：IMU 只解 N300 FDILink；dlog 只认标准 `*.dorec`；无法按图上时段切窗。
+- **改成**：
+  - 新增 `tools/rscap_v2/hi13_imu.py`（HI91：g→m/s²、°/s→rad/s、设备 ms）。
+  - `h32_dlog` 支持 recovered zip（`indices.log` + `data.bin`），ZIP_STORED 成员按文件绝对 offset 直读。
+  - `export_rscap_to_v1.py`：`--imu-kind hi13|n300|auto`、多段 `--imu-rscap`、`--host-start/end` 切窗。
+  - 辅助脚本 `tools/export_usable_20260808_windows.py` 导出优先运动段。
+- **未推送**（按用户要求本地改完即可）。
+
+---
+
 ## 2026-08-05 09:00 (UTC+8)
 
 ### 导出：支持 H32 DLogCapture（MSOP+DIFOP）→ V1

imu_lidar/cli.py

@@ -3,12 +3,82 @@
 from __future__ import annotations
 
 import argparse
+from datetime import datetime
 from pathlib import Path
+from typing import Any
 
 from .contracts import CalibrationMode, CalibrationRequest, CalibrationStatus, SessionInput
+from .phase_a_replay import run_phase_a_replay
 from .pipeline import describe_pipeline, run_calibration
 
 
+def _format_progress_value(value: Any) -> str:
+    if isinstance(value, float):
+        return f"{value:.3f}"
+    if isinstance(value, (list, tuple, set)):
+        return "[" + ",".join(str(item) for item in value) + "]"
+    return str(value)
+
+
+def _print_progress(event: dict[str, Any]) -> None:
+    """Print one compact, immediately flushed progress line."""
+
+    timestamp = datetime.now().strftime("%H:%M:%S")
+    stage_index = event.get("stage_index", "?")
+    stage_total = event.get("stage_total", "?")
+    stage_name = event.get("stage", "unknown")
+    message = event.get("event", "progress")
+    fields = " ".join(
+        f"{key}={_format_progress_value(value)}"
+        for key, value in event.items()
+        if key not in {"stage_index", "stage_total", "stage", "event"}
+        and value is not None
+    )
+    suffix = f" | {fields}" if fields else ""
+    print(
+        f"[{timestamp}] [stage {stage_index}/{stage_total} {stage_name}] {message}{suffix}",
+        flush=True,
+    )
+
+
+def _parse_session_imu_specs(
+    specs: list[str] | None,
+) -> dict[str, Path]:
+    result: dict[str, Path] = {}
+    for spec in specs or []:
+        if "=" not in spec:
+            raise SystemExit(
+                "--session-imu must use SESSION_ID=PATH syntax"
+            )
+        session_id, raw_path = spec.split("=", 1)
+        session_id = session_id.strip()
+        if not session_id or not raw_path.strip():
+            raise SystemExit(
+                "--session-imu must use non-empty SESSION_ID=PATH"
+            )
+        if session_id in result:
+            raise SystemExit(
+                f"duplicate --session-imu for {session_id}"
+            )
+        result[session_id] = Path(raw_path.strip())
+    return result
+
+
+def _print_phase_a_progress(
+    event: str,
+    fields: dict[str, Any],
+) -> None:
+    _print_progress(
+        {
+            "stage_index": "A",
+            "stage_total": "A",
+            "stage": "phase_a_replay",
+            "event": event,
+            **fields,
+        }
+    )
+
+
 def build_parser() -> argparse.ArgumentParser:
     parser = argparse.ArgumentParser(description="LiDAR–IMU extrinsic calibration (V1)")
     subcommands = parser.add_subparsers(dest="command", required=True)
@@ -21,10 +91,28 @@ def build_parser() -> argparse.ArgumentParser:
         default=CalibrationMode.ROTATION_ONLY.value,
     )
 
-    run = subcommands.add_parser("run", help="执行 V1 标定流水线")
-    run.add_argument("--session-id", default="session0")
-    run.add_argument("--imu", required=True, help="IMU CSV/NPZ 路径")
-    run.add_argument("--lidar", required=True, help="LiDAR 会话目录（含 frames_index.csv）")
+    run = subcommands.add_parser(
+        "run",
+        help="执行 V1 标定流水线（可重复 --imu/--lidar/--session-id 做多会话联合）",
+    )
+    run.add_argument(
+        "--session-id",
+        action="append",
+        default=None,
+        help="会话 ID（可重复；与 --imu/--lidar 一一对应）",
+    )
+    run.add_argument(
+        "--imu",
+        action="append",
+        required=True,
+        help="IMU CSV/NPZ 路径（可重复）",
+    )
+    run.add_argument(
+        "--lidar",
+        action="append",
+        required=True,
+        help="LiDAR 会话目录（可重复）",
+    )
     run.add_argument("--vehicle-config", required=True, help="车辆配置 YAML")
     run.add_argument("--output", required=True, help="输出目录")
     run.add_argument(
@@ -34,11 +122,131 @@ def build_parser() -> argparse.ArgumentParser:
     )
     run.add_argument("--max-iterations", type=int, default=2)
     run.add_argument("--time-offset-search-s", type=float, default=1.0)
+    run.add_argument(
+        "--fixed-time-offset-s",
+        type=float,
+        default=None,
+        help="Skip |ω| δt search and use this constant (use 0 after host-UTC bridge)",
+    )
+    run.add_argument(
+        "--session-time-offset-s",
+        action="append",
+        type=float,
+        default=None,
+        help="Per-session fixed time offset; repeat once per --imu/--lidar input",
+    )
+    run.add_argument(
+        "--no-signed-time-refine",
+        action="store_true",
+        help="Disable signed 3-axis δt refine after hand-eye (recommended for host-bridged data)",
+    )
+    run.add_argument(
+        "--max-signed-refine-shift-s",
+        type=float,
+        default=0.05,
+        help="Max |Δδt| accepted by signed refine from the coarse estimate",
+    )
     run.add_argument("--min-pair-rotation-deg", type=float, default=3.0)
     run.add_argument("--min-pair-translation-m", type=float, default=0.3)
+    run.add_argument("--min-registration-fitness", type=float, default=0.5)
+    run.add_argument("--max-imu-gap-s", type=float, default=0.05)
+    run.add_argument("--max-lidar-gap-s", type=float, default=1.0)
+
+    replay = subcommands.add_parser(
+        "phase-a-replay",
+        help="Replay Phase-A from cached motion pairs without rerunning GICP",
+    )
+    replay.add_argument("--motion-pairs", type=Path, required=True)
+    replay.add_argument("--vehicle-config", type=Path, required=True)
+    replay.add_argument("--output", type=Path, required=True)
+    replay.add_argument(
+        "--session-imu",
+        action="append",
+        default=None,
+        metavar="SESSION_ID=PATH",
+        help="Raw IMU mapping used only when cache lacks J_bg/cov",
+    )
+    replay.add_argument(
+        "--exclude-session",
+        action="append",
+        default=None,
+        help="Session ID to exclude; may be repeated",
+    )
+    replay.add_argument(
+        "--strong-rotation-min-deg",
+        type=float,
+        default=1.0,
+    )
+    replay.add_argument(
+        "--decorrelation-block-s",
+        type=float,
+        default=3.0,
+        help="Per-session time-block length used to decorrelate factors",
+    )
+    replay.add_argument(
+        "--max-pairs-per-block",
+        type=int,
+        default=1,
+        help="Maximum factors kept in each decorrelation block",
+    )
+    replay.add_argument(
+        "--bias-prior-sigma-rad-s",
+        type=float,
+        default=0.002,
+    )
+    replay.add_argument(
+        "--yaw-std-max-deg",
+        type=float,
+        default=0.5,
+    )
+    replay.add_argument(
+        "--loo-yaw-range-max-deg",
+        type=float,
+        default=1.0,
+    )
+    replay.add_argument(
+        "--data-prior-difference-max-deg",
+        type=float,
+        default=1.0,
+    )
+    replay.add_argument("--max-nfev", type=int, default=200)
+
     return parser
 
 
+def _build_sessions(args: argparse.Namespace) -> tuple[SessionInput, ...]:
+    imus = [Path(p) for p in args.imu]
+    lidars = [Path(p) for p in args.lidar]
+    if len(imus) != len(lidars):
+        raise SystemExit(f"--imu count ({len(imus)}) must match --lidar count ({len(lidars)})")
+    if args.session_id is None:
+        session_ids = [f"session{i}" for i in range(len(imus))]
+    else:
+        session_ids = list(args.session_id)
+        if len(session_ids) != len(imus):
+            raise SystemExit(
+                f"--session-id count ({len(session_ids)}) must match --imu/--lidar ({len(imus)})"
+            )
+    if args.session_time_offset_s is None:
+        session_offsets: list[float | None] = [None] * len(imus)
+    else:
+        session_offsets = list(args.session_time_offset_s)
+        if len(session_offsets) != len(imus):
+            raise SystemExit(
+                f"--session-time-offset-s count ({len(session_offsets)}) must match "
+                f"--imu/--lidar ({len(imus)})"
+            )
+    return tuple(
+        SessionInput(
+            session_id=sid,
+            imu_source=imu,
+            lidar_source=lidar,
+            fixed_time_offset_s=offset,
+        )
+        for sid, imu, lidar, offset in zip(session_ids, imus, lidars, session_offsets)
+    )
+
+
 def main(argv: list[str] | None = None) -> int:
     parser = build_parser()
     args = parser.parse_args(argv)
@@ -54,28 +262,75 @@ def main(argv: list[str] | None = None) -> int:
             print(f"{index}. {stage.name}: {stage.responsibility}")
         return 0
 
+    if args.command == "phase-a-replay":
+        summary = run_phase_a_replay(
+            motion_pairs_path=args.motion_pairs,
+            vehicle_config_path=args.vehicle_config,
+            output_directory=args.output,
+            imu_paths_by_session=_parse_session_imu_specs(
+                args.session_imu
+            ),
+            excluded_sessions=set(args.exclude_session or []),
+            strong_rotation_min_deg=args.strong_rotation_min_deg,
+            decorrelation_block_s=args.decorrelation_block_s,
+            max_pairs_per_block=args.max_pairs_per_block,
+            bias_prior_sigma_rad_s=args.bias_prior_sigma_rad_s,
+            yaw_std_max_deg=args.yaw_std_max_deg,
+            leave_one_out_yaw_range_max_deg=(
+                args.loo_yaw_range_max_deg
+            ),
+            data_prior_difference_max_deg=(
+                args.data_prior_difference_max_deg
+            ),
+            max_nfev=args.max_nfev,
+            progress_callback=_print_phase_a_progress,
+        )
+        print(f"status: {summary['status']}")
+        print(f"acceptance_checks: {summary['acceptance_checks']}")
+        for name, variant in summary["variants"].items():
+            print(
+                f"{name}: rpy_deg_xyz={variant['rpy_deg_xyz']} "
+                f"RMS={variant['residual_rms_deg']:.6f} "
+                f"P95={variant['residual_p95_deg']:.6f}"
+            )
+        print(
+            "A1 marginalized yaw_std_deg: "
+            f"{summary['marginal_observability_A1']['yaw_std_deg']}"
+        )
+        print(
+            "leave_one_out_yaw_range_deg: "
+            f"{summary['leave_one_out_yaw_range_deg']}"
+        )
+        print(f"report directory: {args.output}")
+        return 0 if (summary["accepted"] or summary.get("partial_accepted")) else 2
+
     if args.command == "run":
+        sessions = _build_sessions(args)
         request = CalibrationRequest(
             vehicle_config=Path(args.vehicle_config),
-            sessions=(
-                SessionInput(
-                    session_id=args.session_id,
-                    imu_source=Path(args.imu),
-                    lidar_source=Path(args.lidar),
-                ),
-            ),
+            sessions=sessions,
             requested_mode=CalibrationMode(args.mode),
             output_directory=Path(args.output),
             max_iterations=args.max_iterations,
             min_pair_rotation_deg=args.min_pair_rotation_deg,
             min_pair_translation_m=args.min_pair_translation_m,
+            min_registration_fitness=args.min_registration_fitness,
+            max_imu_gap_s=args.max_imu_gap_s,
+            max_lidar_gap_s=args.max_lidar_gap_s,
             time_offset_search_s=args.time_offset_search_s,
+            fixed_time_offset_s=args.fixed_time_offset_s,
+            enable_signed_time_refine=not args.no_signed_time_refine,
+            max_signed_refine_shift_s=args.max_signed_refine_shift_s,
         )
-        result = run_calibration(request)
+        result = run_calibration(request, progress_callback=_print_progress)
         print(f"status: {result.status.value}")
         print(f"message: {result.message}")
         if result.time_offset_s is not None:
-            print(f"time_offset_s (t_imu = t_lidar + dt): {result.time_offset_s:.6f}")
+            print(f"time_offset_s (first session; t_imu = t_lidar + dt): {result.time_offset_s:.6f}")
+        joint = (result.details or {}).get("joint") or {}
+        if joint:
+            print(f"merged_pair_count: {joint.get('merged_pair_count')}")
+            print(f"pair_counts_per_session: {joint.get('pair_counts_per_session')}")
         if result.T_IMU_lidar is not None:
             print("T_IMU_lidar:")
             print(result.T_IMU_lidar)

imu_lidar/contracts.py

@@ -25,6 +25,7 @@ class CalibrationStatus(str, Enum):
     NOT_RUN = "not_run"
     BLOCKED = "blocked"
     ROTATION_ONLY_ACCEPTED = "rotation_only_accepted"
+    ROTATION_ONLY_PRIOR_CONSTRAINED = "rotation_only_prior_constrained"
     FULL_SE3_ACCEPTED = "full_se3_accepted"
     FULL_SE3_REJECTED = "full_se3_rejected_due_to_observability"
 
@@ -37,6 +38,9 @@ class SessionInput:
     imu_source: Path
     lidar_source: Path
     board_configuration_id: str | None = None
+    # Optional session-local override. The request-level value remains a
+    # backward-compatible fallback for batches whose timelines are all aligned.
+    fixed_time_offset_s: float | None = None
 
 
 @dataclass(frozen=True)
@@ -50,7 +54,16 @@ class CalibrationRequest:
     max_iterations: int = 2
     min_pair_rotation_deg: float = 3.0
     min_pair_translation_m: float = 0.3
+    min_registration_fitness: float = 0.5
+    max_imu_gap_s: float = 0.05
+    max_lidar_gap_s: float = 1.0
     time_offset_search_s: float = 1.0
+    # If set, skip |ω| search and use this constant (host-UTC-bridged sessions: 0).
+    fixed_time_offset_s: float | None = None
+    # Signed 3-axis refine after hand-eye; disable for already-bridged timelines.
+    enable_signed_time_refine: bool = True
+    # Reject signed refine steps that walk farther than this from the coarse δt.
+    max_signed_refine_shift_s: float = 0.05
 
 
 @dataclass

imu_lidar/finalize.py

@@ -34,6 +34,7 @@ def finalize_result(
     T_IMU_lidar: np.ndarray | None = None,
     time_offset_s: float | None = None,
     output_directory: Path | None = None,
+    motion_pairs_payload: dict[str, Any] | None = None,
 ) -> CalibrationResult:
     """Build the result envelope and optionally write report files."""
 
@@ -72,5 +73,10 @@ def finalize_result(
                 json.dumps({"delta_t_s": time_offset_s, "definition": "t_imu = t_lidar + delta_t"}, indent=2),
                 encoding="utf-8",
             )
+        if motion_pairs_payload is not None:
+            from .motion_pairs_io import save_motion_pairs
+
+            save_motion_pairs(output_directory / "motion_pairs.json", motion_pairs_payload)
+            summary["motion_pairs_file"] = "motion_pairs.json"
         (output_directory / "summary.json").write_text(json.dumps(summary, indent=2), encoding="utf-8")
     return result

imu_lidar/joint_optimizer.py

@@ -2,7 +2,9 @@
 
 from __future__ import annotations
 
-from dataclasses import dataclass
+from collections.abc import Callable, Mapping
+from dataclasses import dataclass, field
+from typing import Any
 
 import numpy as np
 from scipy.optimize import least_squares
@@ -17,10 +19,26 @@ from .imu_preintegration import (
     residual_whiten_matrix,
 )
 from .observability import ObservabilityReport, analyze_observability
+from .phase_a import phase_a_comparison_to_dict, solve_phase_a_comparison
+from .rotation_handeye import select_strong_rotation_pairs
 
 G_NORM = 9.80665
 
 
+@dataclass(frozen=True)
+class PhaseASessionResult:
+    session_id: str
+    pair_count: int
+    gyro_bias0_rad_s: np.ndarray
+    gyro_bias_rad_s: np.ndarray
+    residual_rms_deg: float
+    residual_median_deg: float
+    residual_p95_deg: float
+    outlier_fraction_gt_5deg: float
+    accepted: bool
+    included_in_final: bool
+
+
 @dataclass(frozen=True)
 class JointExtrinsicResult:
     T_IMU_lidar: np.ndarray
@@ -31,6 +49,10 @@ class JointExtrinsicResult:
     gyro_bias_rad_s: np.ndarray | None = None
     accel_bias_m_s2: np.ndarray | None = None
     gravity_m_s2: np.ndarray | None = None
+    gyro_bias_rad_s_per_session: dict[str, np.ndarray] = field(default_factory=dict)
+    phase_a_sessions: tuple[PhaseASessionResult, ...] = ()
+    phase_a_accepted: bool = False
+    phase_a_comparison: dict[str, Any] = field(default_factory=dict)
     notes: tuple[str, ...] = ()
 
 
@@ -152,21 +174,21 @@ def _build_nav_rotations(
     r_x: np.ndarray,
     t_x: np.ndarray,
 ) -> list[np.ndarray]:
-    """Chain IMU orientations in the first-keyframe nav frame using LiDAR+extrinsic."""
+    """Chain IMU orientations; restart at session/gap boundaries (no cross-link)."""
 
+    del id_to_idx
     rotations = [np.eye(3) for _ in keyframe_ids]
     for k in range(len(keyframe_ids) - 1):
         a = keyframe_ids[k]
         b = keyframe_ids[k + 1]
         pair = consecutive_pairs.get((a, b))
         if pair is None:
-            rotations[k + 1] = rotations[k]
+            # Missing link or new session: start a fresh nav chain.
+            rotations[k + 1] = np.eye(3)
             continue
         t_b = np.zeros(3) if pair.t_B_m is None else np.asarray(pair.t_B_m, dtype=float)
         r_meas, _ = _lidar_to_imu_relative(r_x, t_x, pair.R_B, t_b)
         rotations[k + 1] = orthonormalize_rotation(rotations[k] @ r_meas)
-    # Ensure list indexed by id_to_idx
-    del id_to_idx
     return rotations
 
 
@@ -174,10 +196,14 @@ def _solve_phase_c_se3(
     pairs: list[MotionPair],
     r_x: np.ndarray,
     *,
-    gyro_bias0: np.ndarray,
+    gyro_bias_linearization: np.ndarray,
+    gyro_bias_init: np.ndarray,
     gravity_init: np.ndarray,
     sigma_bg_rw: float = 1.0e-5,
     sigma_ba_rw: float = 1.0e-3,
+    t_init: np.ndarray | None = None,
+    t_prior: np.ndarray | None = None,
+    t_prior_sigma_m: np.ndarray | float | None = None,
 ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, float, float, list[str]]:
     """Keyframe IMU factor optimization for full SE(3)."""
 
@@ -185,21 +211,36 @@ def _solve_phase_c_se3(
     usable = [pair for pair in pairs if pair.t_B_m is not None and "delta_v" in pair.metadata]
     if len(usable) < 3:
         notes.append("phase-C skipped: need pairs with full preintegration metadata")
-        return r_x, np.zeros(3), gravity_init, gyro_bias0, np.zeros(3), 1e9, 1e9, notes
+        t0 = np.zeros(3) if t_init is None else np.asarray(t_init, dtype=float).reshape(3)
+        return r_x, t0, gravity_init, gyro_bias_init, np.zeros(3), 1e9, 1e9, notes
 
-    # Unique keyframes sorted by IMU time.
+    # Keyframes: group by session, sort each session by IMU time (no cross-session chain).
     stamp: dict[int, float] = {}
+    kf_session: dict[int, str] = {}
     for pair in usable:
         stamp[pair.i] = float(pair.metadata.get("t_i_imu_s", pair.t_i_s))
         stamp[pair.j] = float(pair.metadata.get("t_j_imu_s", pair.t_j_s))
-    keyframe_ids = sorted(stamp.keys(), key=lambda kid: stamp[kid])
+        kf_session[pair.i] = pair.session_id
+        kf_session[pair.j] = pair.session_id
+    session_ids = sorted(set(kf_session.values()))
+    keyframe_ids: list[int] = []
+    for sid in session_ids:
+        local = [kid for kid, sess in kf_session.items() if sess == sid]
+        local.sort(key=lambda kid: stamp[kid])
+        keyframe_ids.extend(local)
     k_count = len(keyframe_ids)
     id_to_idx = {kid: idx for idx, kid in enumerate(keyframe_ids)}
 
     consecutive_pairs: dict[tuple[int, int], MotionPair] = {}
     for pair in usable:
+        if kf_session.get(pair.i) != kf_session.get(pair.j):
+            continue
         if id_to_idx[pair.j] == id_to_idx[pair.i] + 1:
             consecutive_pairs[(pair.i, pair.j)] = pair
+    notes.append(
+        f"phase-C multi-session graph: sessions={len(session_ids)}, "
+        f"keyframes={k_count}, consecutive_links={len(consecutive_pairs)}"
+    )
 
     g0 = np.asarray(gravity_init, dtype=float).reshape(3)
     if np.linalg.norm(g0) < 1e-6:
@@ -207,16 +248,26 @@ def _solve_phase_c_se3(
     g0 = g0 * (G_NORM / max(np.linalg.norm(g0), 1e-9))
     basis = _gravity_basis(g0)
     ba0 = np.zeros(3)
-    bg0 = np.asarray(gyro_bias0, dtype=float).reshape(3)
+    bg0 = np.asarray(gyro_bias_linearization, dtype=float).reshape(3)
+    bg_init = np.asarray(gyro_bias_init, dtype=float).reshape(3)
 
     # State: dθ(3), t(3), g_xy(2), v(3K), bg(3K), ba(3K)
     n_v = 3 * k_count
     n_b = 3 * k_count
     dim = 3 + 3 + 2 + n_v + n_b + n_b
     x0 = np.zeros(dim)
+    t0 = np.zeros(3) if t_init is None else np.asarray(t_init, dtype=float).reshape(3)
+    x0[3:6] = t0
+    t_prior_vec = None if t_prior is None else np.asarray(t_prior, dtype=float).reshape(3)
+    if t_prior_sigma_m is None:
+        t_sigma = np.array([0.05, 0.05, 0.05], dtype=float)
+    else:
+        t_sigma = np.asarray(t_prior_sigma_m, dtype=float).reshape(-1)
+        if t_sigma.size == 1:
+            t_sigma = np.full(3, float(t_sigma[0]), dtype=float)
     # velocities start at 0; biases at prior
     for idx in range(k_count):
-        x0[8 + n_v + 3 * idx : 8 + n_v + 3 * idx + 3] = bg0
+        x0[8 + n_v + 3 * idx : 8 + n_v + 3 * idx + 3] = bg_init
 
     whitened = []
     for pair in usable:
@@ -269,18 +320,28 @@ def _solve_phase_c_se3(
             w = np.sqrt(_pair_weight(pair))
             out.append(w * (whiten @ err))
 
-        # Bias random-walk between consecutive keyframes.
+        # Bias random-walk between consecutive keyframes (same session only).
         for k in range(k_count - 1):
-            dt = max(stamp[keyframe_ids[k + 1]] - stamp[keyframe_ids[k]], 1e-3)
+            a = keyframe_ids[k]
+            b = keyframe_ids[k + 1]
+            if kf_session.get(a) != kf_session.get(b):
+                continue
+            dt = max(stamp[b] - stamp[a], 1e-3)
             scale_g = 1.0 / (max(sigma_bg_rw, 1e-8) * np.sqrt(dt))
             scale_a = 1.0 / (max(sigma_ba_rw, 1e-8) * np.sqrt(dt))
             out.append(scale_g * (bgs[k + 1] - bgs[k]))
             out.append(scale_a * (bas[k + 1] - bas[k]))
 
-        # Weak priors: first-keyframe biases and translation magnitude.
-        out.append(50.0 * (bgs[0] - bg0))
-        out.append(20.0 * bas[0])
-        out.append(0.2 * t_opt)  # soft |t| prior ~ meters
+        # Weak priors: first keyframe of each session + CAD/installation translation.
+        for sid in session_ids:
+            first = next(kid for kid in keyframe_ids if kf_session[kid] == sid)
+            idx0 = id_to_idx[first]
+            out.append(50.0 * (bgs[idx0] - bg_init))
+            out.append(20.0 * bas[idx0])
+        if t_prior_vec is not None:
+            out.append((t_opt - t_prior_vec) / np.maximum(t_sigma, 1e-3))
+        else:
+            out.append(0.2 * t_opt)  # soft |t|~0 prior when no CAD prior
         return np.concatenate(out)
 
     # Cap evaluations: Phase-C is high-dimensional; synthetic ICP already dominates runtime.
@@ -320,7 +381,197 @@ def _solve_phase_c_se3(
     return r_opt, t_opt, g_opt, bg_mean, ba_mean, rot_rms, trans_rms, notes
 
 
-def solve_joint_extrinsic(
+def _pair_gyro_bias0(pair: MotionPair, fallback: np.ndarray) -> np.ndarray:
+    raw = pair.metadata.get("gyro_bias0_rad_s")
+    if raw is None:
+        return np.asarray(fallback, dtype=float).reshape(3)
+    return np.asarray(raw, dtype=float).reshape(3)
+
+
+def _phase_a_bias_bases(
+    pairs: list[MotionPair],
+    *,
+    gyro_bias_rad_s: np.ndarray | None,
+    gyro_bias_rad_s_by_session: Mapping[str, np.ndarray] | None,
+) -> dict[str, np.ndarray]:
+    session_ids = sorted({pair.session_id for pair in pairs})
+    scalar = None
+    if gyro_bias_rad_s is not None:
+        scalar = np.asarray(gyro_bias_rad_s, dtype=float).reshape(3)
+    supplied = {} if gyro_bias_rad_s_by_session is None else gyro_bias_rad_s_by_session
+    bases: dict[str, np.ndarray] = {}
+    for sid in session_ids:
+        if sid in supplied:
+            bases[sid] = np.asarray(supplied[sid], dtype=float).reshape(3)
+            continue
+        pair = next(
+            (
+                item
+                for item in pairs
+                if item.session_id == sid and "gyro_bias0_rad_s" in item.metadata
+            ),
+            None,
+        )
+        if pair is not None:
+            bases[sid] = np.asarray(pair.metadata["gyro_bias0_rad_s"], dtype=float).reshape(3)
+        elif scalar is not None:
+            bases[sid] = scalar.copy()
+        else:
+            bases[sid] = np.zeros(3)
+    return bases
+
+
+def _rotation_distribution(errs_deg: list[float]) -> tuple[float, float, float, float, bool]:
+    if not errs_deg:
+        return 1e9, 1e9, 1e9, 1.0, False
+    errs = np.asarray(errs_deg, dtype=float)
+    rms = float(np.sqrt(np.mean(errs**2)))
+    median = float(np.median(errs))
+    p95 = float(np.percentile(errs, 95.0))
+    outlier_fraction = float(np.mean(errs > 5.0))
+    accepted = (
+        len(errs) >= 3
+        and rms < 1.5
+        and median < 0.5
+        and p95 < 1.5
+        and outlier_fraction <= 0.005
+    )
+    return rms, median, p95, outlier_fraction, accepted
+
+
+def _solve_phase_a_rotation(
+    pairs: list[MotionPair],
+    r_seed: np.ndarray,
+    *,
+    bias_bases: Mapping[str, np.ndarray],
+    imu: ImuSeries | None,
+    bias_prior_sigma_rad_s: float,
+    preexcluded_session_ids: set[str] | None = None,
+) -> tuple[
+    np.ndarray,
+    dict[str, np.ndarray],
+    tuple[PhaseASessionResult, ...],
+    list[MotionPair],
+    float,
+    bool,
+    list[str],
+]:
+    notes: list[str] = []
+    all_session_ids = sorted({pair.session_id for pair in pairs})
+    prior_w = 1.0 / max(bias_prior_sigma_rad_s, 1e-4)
+
+    def optimize(
+        active_pairs: list[MotionPair],
+        r0: np.ndarray,
+        bias_seed: Mapping[str, np.ndarray],
+    ) -> tuple[np.ndarray, dict[str, np.ndarray]]:
+        session_ids = sorted({pair.session_id for pair in active_pairs})
+        session_index = {sid: index for index, sid in enumerate(session_ids)}
+        whiten = [residual_whiten_matrix(_pair_cov(pair)) for pair in active_pairs]
+        x0 = np.zeros(3 + 3 * len(session_ids))
+        for sid, index in session_index.items():
+            x0[3 + 3 * index : 6 + 3 * index] = np.asarray(bias_seed[sid], dtype=float)
+
+        def residual(vec: np.ndarray) -> np.ndarray:
+            r_opt = orthonormalize_rotation(so3_exp(vec[:3]) @ r0)
+            out: list[np.ndarray] = []
+            for pair, sqrt_info in zip(active_pairs, whiten):
+                index = session_index[pair.session_id]
+                bias = vec[3 + 3 * index : 6 + 3 * index]
+                base = _pair_gyro_bias0(pair, bias_bases[pair.session_id])
+                delta_r = _corrected_delta_r(
+                    pair, bias - base, imu=imu, bias0=base
+                )
+                out.append(
+                    sqrt_info
+                    @ preintegration_rotation_residual(delta_r, r_opt, pair.R_B)
+                )
+            for sid, index in session_index.items():
+                bias = vec[3 + 3 * index : 6 + 3 * index]
+                out.append(prior_w * (bias - bias_bases[sid]))
+            return np.concatenate(out)
+
+        opt = least_squares(residual, x0, loss="huber", f_scale=1.0, max_nfev=200)
+        r_opt = orthonormalize_rotation(so3_exp(opt.x[:3]) @ r0)
+        biases = {
+            sid: opt.x[3 + 3 * index : 6 + 3 * index].copy()
+            for sid, index in session_index.items()
+        }
+        return r_opt, biases
+
+    def summarize(
+        r_opt: np.ndarray,
+        biases: Mapping[str, np.ndarray],
+        included: set[str],
+    ) -> tuple[PhaseASessionResult, ...]:
+        results: list[PhaseASessionResult] = []
+        for sid in all_session_ids:
+            local_pairs = [pair for pair in pairs if pair.session_id == sid]
+            bias = np.asarray(biases.get(sid, bias_bases[sid]), dtype=float).reshape(3)
+            errs: list[float] = []
+            for pair in local_pairs:
+                base = _pair_gyro_bias0(pair, bias_bases[sid])
+                delta_r = _corrected_delta_r(pair, bias - base, imu=imu, bias0=base)
+                err = preintegration_rotation_residual(delta_r, r_opt, pair.R_B)
+                errs.append(float(np.degrees(np.linalg.norm(err))))
+            rms, median, p95, outlier, accepted = _rotation_distribution(errs)
+            results.append(
+                PhaseASessionResult(
+                    session_id=sid,
+                    pair_count=len(local_pairs),
+                    gyro_bias0_rad_s=np.asarray(bias_bases[sid], dtype=float),
+                    gyro_bias_rad_s=bias,
+                    residual_rms_deg=rms,
+                    residual_median_deg=median,
+                    residual_p95_deg=p95,
+                    outlier_fraction_gt_5deg=outlier,
+                    accepted=accepted,
+                    included_in_final=sid in included,
+                )
+            )
+        return tuple(results)
+
+    if not pairs:
+        return r_seed, dict(bias_bases), (), [], 1e9, False, ["no pairs for phase-A"]
+
+    r_first, biases_first = optimize(pairs, r_seed, bias_bases)
+    first = summarize(r_first, biases_first, set(all_session_ids))
+    accepted_ids = {item.session_id for item in first if item.accepted}
+    preexcluded = set() if preexcluded_session_ids is None else set(preexcluded_session_ids)
+    accepted_ids -= preexcluded
+    active_ids = set(all_session_ids)
+    r_final = r_first
+    biases_final = dict(biases_first)
+    if preexcluded and not accepted_ids:
+        active_ids = set()
+        notes.append(f"phase-A pre-gate excluded all sessions: {sorted(preexcluded)}")
+    elif accepted_ids and accepted_ids != active_ids:
+        active_ids = accepted_ids
+        active_pairs = [pair for pair in pairs if pair.session_id in active_ids]
+        r_final, active_biases = optimize(active_pairs, r_first, biases_first)
+        biases_final.update(active_biases)
+        excluded = sorted(set(all_session_ids) - active_ids)
+        notes.append(f"phase-A excluded sessions after local/pre residual gate: {excluded}")
+    active_pairs = [pair for pair in pairs if pair.session_id in active_ids]
+    final = summarize(r_final, biases_final, active_ids)
+    active_results = [item for item in final if item.included_in_final]
+    global_errs: list[float] = []
+    for pair in active_pairs:
+        bias = biases_final[pair.session_id]
+        base = _pair_gyro_bias0(pair, bias_bases[pair.session_id])
+        delta_r = _corrected_delta_r(pair, bias - base, imu=imu, bias0=base)
+        err = preintegration_rotation_residual(delta_r, r_final, pair.R_B)
+        global_errs.append(float(np.degrees(np.linalg.norm(err))))
+    rot_rms, _, _, _, global_ok = _rotation_distribution(global_errs)
+    accepted = bool(active_results and global_ok and all(item.accepted for item in active_results))
+    notes.append(
+        f"phase-A session-local bias refine: sessions={len(active_ids)}/{len(all_session_ids)}, "
+        f"pairs={len(active_pairs)}, rms={rot_rms:.3f} deg"
+    )
+    return r_final, biases_final, final, active_pairs, rot_rms, accepted, notes
+
+
+def _solve_joint_extrinsic_legacy(
     pairs: list[MotionPair] | tuple[MotionPair, ...],
     r_x: np.ndarray,
     *,
@@ -329,8 +580,13 @@ def solve_joint_extrinsic(
     delta_t_s: float = 0.0,
     gyro_bias_rad_s: np.ndarray | None = None,
     gravity_init_m_s2: np.ndarray | None = None,
+    gyro_bias_rad_s_by_session: Mapping[str, np.ndarray] | None = None,
+    time_offset_s_by_session: Mapping[str, float] | None = None,
     bias_prior_sigma_rad_s: float = 0.02,
     enable_phase_c: bool | None = None,
+    t_init_m: np.ndarray | None = None,
+    t_prior_m: np.ndarray | None = None,
+    t_prior_sigma_m: np.ndarray | float | None = None,
 ) -> JointExtrinsicResult:
     """Refine extrinsic using Phase-A whitened rotation factors, optional Phase-C SE(3)."""
 
@@ -344,16 +600,17 @@ def solve_joint_extrinsic(
 
     r = orthonormalize_rotation(np.asarray(r_x, dtype=float))
     bias0 = np.zeros(3) if gyro_bias_rad_s is None else np.asarray(gyro_bias_rad_s, dtype=float).reshape(3)
+    t_seed = None if t_init_m is None else np.asarray(t_init_m, dtype=float).reshape(3)
     weights = np.asarray([_pair_weight(pair) for pair in usable], dtype=float)
     whitens = [residual_whiten_matrix(_pair_cov(pair)) for pair in usable]
     prior_w = 1.0 / max(bias_prior_sigma_rad_s, 1e-4)
 
     def rotation_residuals(r_opt: np.ndarray, delta_bias: np.ndarray) -> np.ndarray:
         residuals = []
-        for pair, weight, whiten in zip(usable, weights, whitens):
+        for pair, whiten in zip(usable, whitens):
             delta_r = _corrected_delta_r(pair, delta_bias, imu=imu, bias0=bias0)
             err = preintegration_rotation_residual(delta_r, r_opt, pair.R_B)
-            residuals.append(np.sqrt(weight) * (whiten @ err))
+            residuals.append(whiten @ err)
         residuals.append(prior_w * delta_bias)
         return np.concatenate(residuals) if residuals else np.zeros(0)
 
@@ -366,16 +623,16 @@ def solve_joint_extrinsic(
             residual_rot_bias,
             np.zeros(6),
             loss="huber",
-            f_scale=np.deg2rad(1.0),
+            f_scale=1.0,
             max_nfev=200,
         )
         r = orthonormalize_rotation(so3_exp(opt.x[:3]) @ r)
         delta_bias = opt.x[3:]
         bias_out = bias0 + delta_bias
         notes.append(
-            "phase-A joint refine (Σ-whitened + J_bg): "
+            "phase-A joint refine (single Σ whitening + J_bg): "
             f"|δb|={float(np.linalg.norm(delta_bias)):.3e} rad/s, "
-            f"weighted pairs={len(usable)}"
+            f"pairs={len(usable)}"
         )
     else:
         bias_out = bias0
@@ -389,7 +646,7 @@ def solve_joint_extrinsic(
         rot_errs.append(np.degrees(np.linalg.norm(err)))
     rot_rms = float(np.sqrt(np.mean(np.square(rot_errs)))) if rot_errs else 1e9
 
-    t = np.zeros(3)
+    t = np.zeros(3) if t_seed is None else t_seed.copy()
     translation_accepted = False
     trans_rms = 1e9
     gravity_out: np.ndarray | None = None
@@ -400,6 +657,12 @@ def solve_joint_extrinsic(
     else:
         gravity_init = np.asarray(gravity_init_m_s2, dtype=float).reshape(3)
 
+    if t_prior_m is not None:
+        notes.append(
+            "using CAD/installation translation prior "
+            f"t={np.asarray(t_prior_m, dtype=float).reshape(3).tolist()}"
+        )
+
     if (
         enable_phase_c
         and not force_rotation_only
@@ -410,14 +673,26 @@ def solve_joint_extrinsic(
         r, t, gravity_out, bias_out, accel_bias_out, rot_rms, trans_rms, c_notes = _solve_phase_c_se3(
             usable,
             r,
-            gyro_bias0=bias_out,
+            gyro_bias_linearization=bias0,
+            gyro_bias_init=bias_out,
             gravity_init=gravity_init,
+            t_init=t_seed if t_seed is not None else t_prior_m,
+            t_prior=t_prior_m,
+            t_prior_sigma_m=t_prior_sigma_m,
         )
         notes.extend(c_notes)
         translation_accepted = bool(trans_rms < 0.75 and np.linalg.norm(t) > 1e-4)
         if not translation_accepted:
-            notes.append("phase-C translation residual/gate failed; keeping translation at zero")
-            t = np.zeros(3)
+            # Prefer CAD prior over silent zero when motion SE3 is rejected.
+            if t_prior_m is not None:
+                t = np.asarray(t_prior_m, dtype=float).reshape(3)
+                notes.append(
+                    "phase-C translation residual/gate failed; CAD translation is reported "
+                    "as a prior only and is not accepted as calibration"
+                )
+            else:
+                notes.append("phase-C translation residual/gate failed; keeping translation at zero")
+                t = np.zeros(3)
     elif (
         not force_rotation_only
         and observability.translation_observable
@@ -437,9 +712,19 @@ def solve_joint_extrinsic(
                 pred = (pair.R_A - np.eye(3)) @ t_opt
                 meas = r_opt @ np.asarray(pair.t_B_m, dtype=float)
                 residuals.append(np.sqrt(weight) * (pred - meas))
+            if t_prior_m is not None:
+                sigma = np.asarray(t_prior_sigma_m if t_prior_sigma_m is not None else 0.05, dtype=float)
+                if sigma.size == 1:
+                    sigma = np.full(3, float(sigma), dtype=float)
+                residuals.append((t_opt - np.asarray(t_prior_m, dtype=float).reshape(3)) / np.maximum(sigma, 1e-3))
             return np.concatenate(residuals)
 
-        opt_t = least_squares(residual_se3, np.zeros(6), loss="huber", f_scale=0.05, max_nfev=200)
+        x_se3 = np.zeros(6)
+        if t_seed is not None:
+            x_se3[3:] = t_seed
+        elif t_prior_m is not None:
+            x_se3[3:] = np.asarray(t_prior_m, dtype=float).reshape(3)
+        opt_t = least_squares(residual_se3, x_se3, loss="huber", f_scale=0.05, max_nfev=200)
         r = orthonormalize_rotation(so3_exp(opt_t.x[:3]) @ r)
         t = opt_t.x[3:]
         rot_errs = []
@@ -457,6 +742,13 @@ def solve_joint_extrinsic(
         if not translation_accepted:
             notes.append("translation residual too large; keeping translation at zero")
             t = np.zeros(3)
+    elif not force_rotation_only and t_prior_m is not None:
+        t = np.asarray(t_prior_m, dtype=float).reshape(3)
+        translation_accepted = False
+        notes.append(
+            "SE3 motion solve gated off; CAD translation is reported as a prior only "
+            "and is not accepted as calibration"
+        )
     else:
         notes.append("rotation-only extrinsic returned (phase-A; phase-C SE3 gated off)")
 
@@ -464,10 +756,222 @@ def solve_joint_extrinsic(
         T_IMU_lidar=make_transform(t, r),
         translation_accepted=bool(translation_accepted and np.linalg.norm(t) > 0),
         residual_rms_rot_deg=rot_rms,
-        residual_rms_trans_m=0.0 if not translation_accepted else trans_rms,
+        residual_rms_trans_m=trans_rms,
         observability=observability,
         gyro_bias_rad_s=np.asarray(bias_out, dtype=float),
         accel_bias_m_s2=None if accel_bias_out is None else np.asarray(accel_bias_out, dtype=float),
         gravity_m_s2=None if gravity_out is None else np.asarray(gravity_out, dtype=float),
         notes=tuple(notes),
     )
+
+
+def solve_joint_extrinsic(
+    pairs: list[MotionPair] | tuple[MotionPair, ...],
+    r_x: np.ndarray,
+    *,
+    force_rotation_only: bool = False,
+    imu: ImuSeries | None = None,
+    delta_t_s: float = 0.0,
+    gyro_bias_rad_s: np.ndarray | None = None,
+    gyro_bias_rad_s_by_session: Mapping[str, np.ndarray] | None = None,
+    time_offset_s_by_session: Mapping[str, float] | None = None,
+    preexcluded_session_ids: set[str] | None = None,
+    gravity_init_m_s2: np.ndarray | None = None,
+    bias_prior_sigma_rad_s: float = 0.002,
+    rotation_prior: np.ndarray | None = None,
+    rotation_prior_sigma_deg: float = 15.0,
+    phase_a_yaw_std_max_deg: float = 0.5,
+    phase_a_loo_yaw_range_max_deg: float = 1.0,
+    phase_a_data_prior_difference_max_deg: float = 1.0,
+    run_phase_a_leave_one_out: bool = True,
+    phase_a_progress_callback: (
+        Callable[[str, dict[str, Any]], None] | None
+    ) = None,
+    enable_phase_c: bool | None = None,
+    t_init_m: np.ndarray | None = None,
+    t_prior_m: np.ndarray | None = None,
+    t_prior_sigma_m: np.ndarray | float | None = None,
+) -> JointExtrinsicResult:
+    """Run the corrected session-aware Phase-A and gate unfinished SE(3) stages."""
+
+    del gravity_init_m_s2, t_init_m, t_prior_sigma_m, imu, r_x
+    usable_input = [pair for pair in pairs if pair.t_B_m is not None]
+    bias_bases = _phase_a_bias_bases(
+        usable_input,
+        gyro_bias_rad_s=gyro_bias_rad_s,
+        gyro_bias_rad_s_by_session=gyro_bias_rad_s_by_session,
+    )
+    comparison = solve_phase_a_comparison(
+        usable_input,
+        gyro_bias_rad_s_by_session=bias_bases,
+        rotation_prior=rotation_prior,
+        rotation_prior_sigma_deg=rotation_prior_sigma_deg,
+        preexcluded_session_ids=preexcluded_session_ids,
+        bias_prior_sigma_rad_s=bias_prior_sigma_rad_s,
+        yaw_std_max_deg=phase_a_yaw_std_max_deg,
+        leave_one_out_yaw_range_max_deg=(
+            phase_a_loo_yaw_range_max_deg
+        ),
+        data_prior_difference_max_deg=(
+            phase_a_data_prior_difference_max_deg
+        ),
+        run_leave_one_out=run_phase_a_leave_one_out,
+        progress_callback=phase_a_progress_callback,
+    )
+    primary = comparison.session_bg_data_only
+    r = primary.R_IMU_lidar
+    biases = primary.gyro_bias_rad_s_per_session
+    rot_rms = primary.residual_rms_deg
+    phase_a_accepted = comparison.accepted
+    notes = list(comparison.notes)
+    notes.append(
+        "phase-A primary=A1_session_bg_data_only; "
+        f"A0 RPY={comparison.fixed_bg_data_only.rpy_deg_xyz.tolist()}, "
+        f"A1 RPY={primary.rpy_deg_xyz.tolist()}, "
+        "A2 RPY="
+        f"{comparison.session_bg_with_rotation_prior.rpy_deg_xyz.tolist()}"
+    )
+    notes.append(
+        f"phase-A marginal yaw_std={comparison.marginal_observability.yaw_std_deg:.3f} deg, "
+        f"LOO yaw range={comparison.leave_one_out_yaw_range_deg:.3f} deg"
+    )
+
+    session_results_list: list[PhaseASessionResult] = [
+        PhaseASessionResult(
+            session_id=item.session_id,
+            pair_count=item.pair_count,
+            gyro_bias0_rad_s=item.gyro_bias0_rad_s,
+            gyro_bias_rad_s=item.gyro_bias_rad_s,
+            residual_rms_deg=item.residual_rms_deg,
+            residual_median_deg=item.residual_median_deg,
+            residual_p95_deg=item.residual_p95_deg,
+            outlier_fraction_gt_5deg=item.outlier_fraction_gt_5deg,
+            accepted=item.accepted,
+            included_in_final=True,
+        )
+        for item in primary.sessions
+    ]
+    preexcluded = (
+        set()
+        if preexcluded_session_ids is None
+        else set(preexcluded_session_ids)
+    )
+    strong_all = select_strong_rotation_pairs(usable_input)
+    for session_id in sorted(preexcluded):
+        local_pairs = [
+            pair for pair in strong_all if pair.session_id == session_id
+        ]
+        errors = [
+            float(
+                np.degrees(
+                    np.linalg.norm(
+                        preintegration_rotation_residual(
+                            pair.R_A, r, pair.R_B
+                        )
+                    )
+                )
+            )
+            for pair in local_pairs
+        ]
+        rms, median, p95, outlier, accepted = _rotation_distribution(
+            errors
+        )
+        base = np.asarray(
+            bias_bases.get(session_id, np.zeros(3)), dtype=float
+        ).reshape(3)
+        session_results_list.append(
+            PhaseASessionResult(
+                session_id=session_id,
+                pair_count=len(local_pairs),
+                gyro_bias0_rad_s=base,
+                gyro_bias_rad_s=base,
+                residual_rms_deg=rms,
+                residual_median_deg=median,
+                residual_p95_deg=p95,
+                outlier_fraction_gt_5deg=outlier,
+                accepted=accepted,
+                included_in_final=False,
+            )
+        )
+    session_results = tuple(
+        sorted(session_results_list, key=lambda item: item.session_id)
+    )
+    usable = [
+        pair
+        for pair in strong_all
+        if pair.session_id not in preexcluded
+    ]
+    base_observability = analyze_observability(usable, r)
+    marginal = comparison.marginal_observability
+    observability = ObservabilityReport(
+        rotation_observable=bool(
+            marginal.rank == 3
+            and marginal.yaw_std_deg <= phase_a_yaw_std_max_deg
+        ),
+        translation_observable=base_observability.translation_observable,
+        condition_rotation=marginal.condition,
+        condition_translation=base_observability.condition_translation,
+        notes=tuple(
+            list(marginal.notes)
+            + list(base_observability.notes)
+        ),
+    )
+    notes.extend(observability.notes)
+    if time_offset_s_by_session is None:
+        notes.append(
+            f"legacy scalar time offset fixed during pair construction: {float(delta_t_s):.6f}s"
+        )
+    else:
+        fixed_offsets = {
+            str(sid): float(value) for sid, value in time_offset_s_by_session.items()
+        }
+        notes.append(
+            f"time offsets fixed during pair construction (not optimized): {fixed_offsets}"
+        )
+
+    for item in session_results:
+        notes.append(
+            f"phase-A session {item.session_id}: included={item.included_in_final}, "
+            f"pairs={item.pair_count}, rms={item.residual_rms_deg:.3f} deg, "
+            f"p95={item.residual_p95_deg:.3f} deg, "
+            f"|bias-bias0|={float(np.linalg.norm(item.gyro_bias_rad_s - item.gyro_bias0_rad_s)):.3e}"
+        )
+
+    phase_c_requested = (not force_rotation_only) if enable_phase_c is None else bool(enable_phase_c)
+    t = np.zeros(3)
+    if not force_rotation_only:
+        if phase_c_requested:
+            notes.append(
+                "phase-B/C gated off: session-aware translation/gravity/navigation "
+                "states are not implemented yet"
+            )
+        else:
+            notes.append("phase-C disabled; translation is not accepted")
+        if t_prior_m is not None:
+            t = np.asarray(t_prior_m, dtype=float).reshape(3)
+            notes.append(
+                "CAD translation is reported as a prior only and is not accepted as calibration"
+            )
+    else:
+        notes.append("rotation-only extrinsic returned after corrected phase-A")
+
+    single_bias = None
+    if len(biases) == 1:
+        single_bias = np.asarray(next(iter(biases.values())), dtype=float)
+    return JointExtrinsicResult(
+        T_IMU_lidar=make_transform(t, r),
+        translation_accepted=False,
+        residual_rms_rot_deg=rot_rms,
+        residual_rms_trans_m=1e9,
+        observability=observability,
+        gyro_bias_rad_s=single_bias,
+        accel_bias_m_s2=None,
+        gravity_m_s2=None,
+        gyro_bias_rad_s_per_session={
+            sid: np.asarray(value, dtype=float) for sid, value in biases.items()
+        },
+        phase_a_sessions=session_results,
+        phase_a_accepted=phase_a_accepted,
+        phase_a_comparison=phase_a_comparison_to_dict(comparison),
+        notes=tuple(notes),
+    )

imu_lidar/keyframes.py

@@ -21,6 +21,7 @@ def build_keyframes(
     *,
     min_translation_m: float = 0.3,
     min_rotation_deg: float = 3.0,
+    min_registration_fitness: float = 0.5,
     max_frame_gap: int = 8,
 ) -> KeyframeSet:
     """Select keyframes with enough relative motion for hand-eye pairs."""
@@ -36,7 +37,7 @@ def build_keyframes(
             last = index
             continue
         result = register_lidar_pair(frames[index].points_xyz, frames[last].points_xyz)
-        if not result.ok:
+        if not result.ok or result.fitness < min_registration_fitness:
             continue
         if result.translation_m >= min_translation_m or result.rotation_deg >= min_rotation_deg:
             selected.append(index)

imu_lidar/lidar_io.py

@@ -19,10 +19,7 @@ import numpy as np
 from .contracts import LidarFrame
 
 
-def load_lidar_frames(path: Path | str) -> list[LidarFrame]:
-    """Load all LiDAR frames listed by ``frames_index.csv`` under ``path``."""
-
-    root = Path(path)
+def _read_frames_index(root: Path) -> tuple[np.ndarray, str]:
     index_path = root / "frames_index.csv"
     if not index_path.exists():
         raise FileNotFoundError(f"missing frames_index.csv under {root}")
@@ -38,6 +35,54 @@ def load_lidar_frames(path: Path | str) -> list[LidarFrame]:
         raise ValueError(
             f"frames_index.csv must contain frame_id,{file_key}/filename,t_start,t_end; got {sorted(names)}"
         )
+    return rows, file_key
+
+
+def list_lidar_frame_entries(path: Path | str) -> list[tuple[str, float, float, Path]]:
+    """Return ``(frame_id, t_start, t_end, npz_path)`` sorted by mid time (same as ``load_lidar_frames``)."""
+
+    root = Path(path)
+    rows, file_key = _read_frames_index(root)
+    entries: list[tuple[str, float, float, Path]] = []
+    for row in rows:
+        t0 = float(row["t_start"])
+        t1 = float(row["t_end"])
+        entries.append((str(row["frame_id"]), t0, t1, root / str(row[file_key])))
+    entries.sort(key=lambda item: 0.5 * (item[1] + item[2]))
+    return entries
+
+
+def load_lidar_frame_at(root: Path | str, index: int) -> LidarFrame:
+    """Load one frame by index in mid-time-sorted order (matches motion-pair ``i``/``j``)."""
+
+    entries = list_lidar_frame_entries(root)
+    if index < 0 or index >= len(entries):
+        raise IndexError(f"frame index {index} outside [0, {len(entries) - 1}] for {root}")
+    frame_id, t0, t1, npz_path = entries[index]
+    with np.load(npz_path) as payload:
+        if "points" not in payload.files:
+            raise ValueError(f"{npz_path} must contain array 'points'")
+        points = np.asarray(payload["points"], dtype=float)
+    if points.ndim != 2 or points.shape[1] < 3:
+        raise ValueError(f"{npz_path}: points must have shape (N, 3[+])")
+    return LidarFrame(
+        frame_id=frame_id,
+        t_start_s=t0,
+        t_end_s=t1,
+        points_xyz=points[:, :3],
+        path=npz_path,
+    )
+
+
+def lidar_frame_count(path: Path | str) -> int:
+    return len(list_lidar_frame_entries(path))
+
+
+def load_lidar_frames(path: Path | str) -> list[LidarFrame]:
+    """Load all LiDAR frames listed by ``frames_index.csv`` under ``path``."""
+
+    root = Path(path)
+    rows, file_key = _read_frames_index(root)
 
     frames: list[LidarFrame] = []
     for row in rows:

imu_lidar/motion_pairs.py

@@ -2,7 +2,10 @@
 
 from __future__ import annotations
 
+from collections.abc import Callable
 from dataclasses import dataclass
+from time import perf_counter
+from typing import Any
 
 import numpy as np
 
@@ -30,7 +33,12 @@ def build_motion_pairs(
     acc_bias_m_s2: np.ndarray | None = None,
     min_rotation_deg: float = 3.0,
     min_translation_m: float = 0.3,
+    min_registration_fitness: float = 0.5,
+    max_imu_gap_s: float = 0.05,
+    max_lidar_gap_s: float = 1.0,
+    all_frame_times_s: np.ndarray | None = None,
     max_index_span: int = 4,
+    progress_callback: Callable[[dict[str, Any]], None] | None = None,
 ) -> MotionPairSet:
     """Create A/B motion pairs between nearby keyframes.
 
@@ -40,21 +48,74 @@ def build_motion_pairs(
 
     notes: list[str] = []
     pairs: list[MotionPair] = []
+    rejected_fitness = 0
+    rejected_imu_gap = 0
+    rejected_lidar_gap = 0
+    frame_times = (
+        None
+        if all_frame_times_s is None
+        else np.asarray(all_frame_times_s, dtype=float).reshape(-1)
+    )
     bias_g = np.zeros(3) if gyro_bias_rad_s is None else np.asarray(gyro_bias_rad_s, dtype=float)
     bias_a = np.zeros(3) if acc_bias_m_s2 is None else np.asarray(acc_bias_m_s2, dtype=float)
     n = len(keyframes)
     if n < 2:
         return MotionPairSet((), ("need at least two keyframes",))
 
+    total_candidates = sum(max(n - span, 0) for span in range(1, max_index_span + 1))
+    processed_candidates = 0
+    started_at = perf_counter()
+    last_progress_at = started_at
+
+    def report_progress(*, event: str, span: int, force: bool = False) -> None:
+        nonlocal last_progress_at
+        if progress_callback is None:
+            return
+        now = perf_counter()
+        if not force and processed_candidates > 1 and now - last_progress_at < 10.0:
+            return
+        last_progress_at = now
+        progress_callback(
+            {
+                "event": event,
+                "processed_candidates": processed_candidates,
+                "total_candidates": total_candidates,
+                "progress_pct": 100.0 * processed_candidates / max(total_candidates, 1),
+                "current_span": span,
+                "max_span": max_index_span,
+                "accepted_pairs": len(pairs),
+                "rejected_fitness": rejected_fitness,
+                "rejected_imu_gap": rejected_imu_gap,
+                "rejected_lidar_gap": rejected_lidar_gap,
+                "elapsed_s": now - started_at,
+            }
+        )
+
+    report_progress(event="start", span=1, force=True)
+
     for span in range(1, max_index_span + 1):
         for start in range(0, n - span):
+            processed_candidates += 1
+            report_progress(event="running", span=span)
             i = start
             j = start + span
             frame_i = keyframes[i]
             frame_j = keyframes[j]
+            source_i = int(keyframe_indices[i])
+            source_j = int(keyframe_indices[j])
+            if frame_times is not None:
+                lo = min(source_i, source_j)
+                hi = max(source_i, source_j)
+                local_times = frame_times[lo : hi + 1]
+                if local_times.size >= 2 and np.any(np.diff(local_times) > max_lidar_gap_s):
+                    rejected_lidar_gap += 1
+                    continue
             reg = register_lidar_pair(frame_j.points_xyz, frame_i.points_xyz)
             if not reg.ok:
                 continue
+            if reg.fitness < min_registration_fitness:
+                rejected_fitness += 1
+                continue
             if reg.rotation_deg < min_rotation_deg and reg.translation_m < min_translation_m:
                 continue
 
@@ -64,6 +125,16 @@ def build_motion_pairs(
                 continue
             if t_i_imu < imu.t_s[0] or t_j_imu > imu.t_s[-1]:
                 continue
+            imu_lo = max(int(np.searchsorted(imu.t_s, t_i_imu, side="right")) - 1, 0)
+            imu_hi = min(
+                int(np.searchsorted(imu.t_s, t_j_imu, side="left")) + 1,
+                imu.t_s.size,
+            )
+            if imu_hi - imu_lo >= 2 and np.any(
+                np.diff(imu.t_s[imu_lo:imu_hi]) > max_imu_gap_s
+            ):
+                rejected_imu_gap += 1
+                continue
 
             preint = preintegrate_imu(
                 imu.t_s,
@@ -111,14 +182,27 @@ def build_motion_pairs(
                         "delta_p": preint.delta_p.tolist(),
                         "t_i_imu_s": t_i_imu,
                         "t_j_imu_s": t_j_imu,
+                        "gyro_bias0_rad_s": bias_g.tolist(),
+                        "accel_bias0_m_s2": bias_a.tolist(),
+                        "time_offset_s": float(delta_t_s),
+                        "keyframe_span": int(span),
+                        "is_consecutive": bool(span == 1),
                         "modeling": "imu_preintegration_factor_phase_c",
                     },
                 )
             )
 
+    report_progress(event="complete", span=max_index_span, force=True)
+
     notes.append(
         f"built {len(pairs)} motion pairs (Phase-C preintegration: ΔR/Δv/Δp, Σ9, J_bg/J_ba)"
     )
+    notes.append(
+        "quality rejects: "
+        f"fitness<{min_registration_fitness:.2f}: {rejected_fitness}, "
+        f"IMU gap>{max_imu_gap_s:.3f}s: {rejected_imu_gap}, "
+        f"LiDAR gap>{max_lidar_gap_s:.3f}s: {rejected_lidar_gap}"
+    )
     return MotionPairSet(pairs=tuple(pairs), notes=tuple(notes))
 
 

imu_lidar/motion_pairs_io.py

@@ -0,0 +1,154 @@
+"""Serialize / deserialize motion pairs for fast visualization."""
+
+from __future__ import annotations
+
+import json
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+
+from .contracts import MotionPair
+
+SCHEMA_VERSION = 2
+
+# Keep visualization fields plus the compact 3x3 rotation metadata needed to
+# rerun Phase-A without repeating LiDAR registration. Full 9x9 Phase-C matrices
+# remain excluded from this cache.
+_METADATA_KEEP = frozenset(
+    {
+        "backend",
+        "rotation_deg_A",
+        "rotation_deg_B",
+        "translation_m_B",
+        "weight",
+        "duration_s",
+        "mean_gyro_norm",
+        "preint_sigma_rad",
+        "cov",
+        "J_bg",
+        "phase_a_metadata_rehydrated",
+        "rehydrated_R_A_error_deg",
+        "t_i_imu_s",
+        "t_j_imu_s",
+        "gyro_bias0_rad_s",
+        "accel_bias0_m_s2",
+        "time_offset_s",
+        "keyframe_span",
+        "is_consecutive",
+        "modeling",
+    }
+)
+
+
+def _to_list(value: Any) -> Any:
+    if isinstance(value, np.ndarray):
+        return value.tolist()
+    if isinstance(value, (np.floating, np.integer, np.bool_)):
+        return value.item()
+    return value
+
+
+def pair_to_dict(pair: MotionPair) -> dict[str, Any]:
+    meta = {
+        str(k): _to_list(v)
+        for k, v in (pair.metadata or {}).items()
+        if str(k) in _METADATA_KEEP
+    }
+    return {
+        "session_id": pair.session_id,
+        "i": int(pair.i),
+        "j": int(pair.j),
+        "t_i_s": float(pair.t_i_s),
+        "t_j_s": float(pair.t_j_s),
+        "R_A": np.asarray(pair.R_A, dtype=float).reshape(3, 3).tolist(),
+        "R_B": np.asarray(pair.R_B, dtype=float).reshape(3, 3).tolist(),
+        "t_A_m": None if pair.t_A_m is None else np.asarray(pair.t_A_m, dtype=float).reshape(3).tolist(),
+        "t_B_m": None if pair.t_B_m is None else np.asarray(pair.t_B_m, dtype=float).reshape(3).tolist(),
+        "fitness": float(pair.fitness),
+        "metadata": meta,
+    }
+
+
+def pair_from_dict(payload: dict[str, Any]) -> MotionPair:
+    t_a = payload.get("t_A_m")
+    t_b = payload.get("t_B_m")
+    return MotionPair(
+        session_id=str(payload.get("session_id", "")),
+        i=int(payload["i"]),
+        j=int(payload["j"]),
+        t_i_s=float(payload["t_i_s"]),
+        t_j_s=float(payload["t_j_s"]),
+        R_A=np.asarray(payload["R_A"], dtype=float).reshape(3, 3),
+        R_B=np.asarray(payload["R_B"], dtype=float).reshape(3, 3),
+        t_A_m=None if t_a is None else np.asarray(t_a, dtype=float).reshape(3),
+        t_B_m=None if t_b is None else np.asarray(t_b, dtype=float).reshape(3),
+        fitness=float(payload.get("fitness", 0.0)),
+        metadata=dict(payload.get("metadata") or {}),
+    )
+
+
+def build_motion_pairs_payload(
+    *,
+    prepared_sessions: list[dict[str, Any]],
+) -> dict[str, Any]:
+    """Build a JSON-serializable cache from pipeline ``prepared`` session dicts."""
+
+    sessions_out: list[dict[str, Any]] = []
+    for prep in prepared_sessions:
+        pairs = prep.get("pairs") or ()
+        sessions_out.append(
+            {
+                "session_id": prep.get("session_id"),
+                "delta_t_s": float(prep.get("time_offset_s", 0.0)),
+                "gyro_bias_rad_s": np.asarray(prep.get("gyro_bias_rad_s", np.zeros(3)), dtype=float)
+                .reshape(3)
+                .tolist(),
+                "pair_count": len(pairs),
+                "pairs": [pair_to_dict(pair) for pair in pairs],
+            }
+        )
+    return {
+        "schema_version": SCHEMA_VERSION,
+        "sessions": sessions_out,
+        "note": "Cached motion pairs for visualization; A=IMU preintegration, B=LiDAR registration",
+    }
+
+
+def save_motion_pairs(path: Path | str, payload: dict[str, Any]) -> Path:
+    destination = Path(path)
+    destination.parent.mkdir(parents=True, exist_ok=True)
+    destination.write_text(json.dumps(payload, indent=2), encoding="utf-8")
+    return destination
+
+
+def load_motion_pairs(path: Path | str) -> dict[str, Any]:
+    payload = json.loads(Path(path).read_text(encoding="utf-8"))
+    version = int(payload.get("schema_version", 0))
+    if version not in {1, SCHEMA_VERSION}:
+        raise ValueError(
+            f"unsupported motion_pairs schema_version={version}; "
+            f"expected 1 or {SCHEMA_VERSION}"
+        )
+    return payload
+
+
+def pairs_for_session(payload: dict[str, Any], session_id: str | None = None) -> list[MotionPair]:
+    sessions = payload.get("sessions") or []
+    if not sessions:
+        return []
+    if session_id is None:
+        chosen = sessions[0]
+    else:
+        chosen = next((s for s in sessions if s.get("session_id") == session_id), None)
+        if chosen is None:
+            raise KeyError(f"session_id {session_id!r} not found in motion_pairs cache")
+    return [pair_from_dict(item) for item in chosen.get("pairs") or []]
+
+
+def resolve_motion_pairs_path(summary_path: Path | str) -> Path | None:
+    """Return ``motion_pairs.json`` next to a summary if it exists."""
+
+    summary = Path(summary_path)
+    candidate = summary.parent / "motion_pairs.json"
+    return candidate if candidate.is_file() else None

imu_lidar/observability.py

@@ -57,43 +57,53 @@ def analyze_observability(
     if j_r.size == 0:
         return ObservabilityReport(False, False, 1e9, 1e9, ("empty rotation jacobian",))
 
-    # Normalize columns.
-    col_norm = np.linalg.norm(j_r, axis=0) + 1e-12
-    j_r_n = j_r / col_norm
-    singular = np.linalg.svd(j_r_n, compute_uv=False)
+    singular = np.linalg.svd(j_r, compute_uv=False)
     cond_r = float(singular[0] / max(singular[-1], 1e-12))
-    rotation_ok = cond_r < condition_threshold and singular[-1] > 1e-3
+    rotation_information = float(singular[-1] / np.sqrt(max(len(usable), 1)))
+    rotation_ok = (
+        cond_r < condition_threshold
+        and rotation_information > 1e-3
+        and singular[-1] > 1e-6
+    )
 
-    # Translation observability proxy: diversity of rotation axes and presence of translation in B.
-    axes = []
-    translations = []
-    for pair in usable:
-        axis = so3_log(pair.R_B)
-        n = np.linalg.norm(axis)
-        if n > 1e-8:
-            axes.append(axis / n)
-        if pair.t_B_m is not None:
-            translations.append(pair.t_B_m)
-    axis_rank = 0
-    if axes:
-        axis_mat = np.asarray(axes, dtype=float)
-        axis_rank = int(np.linalg.matrix_rank(axis_mat, tol=0.1))
-    trans_span = 0.0
-    if translations:
-        tmat = np.asarray(translations, dtype=float)
-        trans_span = float(np.linalg.norm(np.std(tmat, axis=0)))
-    # For planar yaw-mostly motion, translation z is typically weak.
-    translation_ok = axis_rank >= 2 and trans_span > 0.2 and len(translations) >= 5
-    cond_t = 1e9 if not translation_ok else float(max(3, 10 - axis_rank * 2) * (0.5 / max(trans_span, 1e-3)))
+    # Translation lever arm is observable through stacked (R_A - I). Pure
+    # planar yaw leaves its vertical column in the nullspace and must fail.
+    translation_rows = [
+        np.asarray(pair.R_A, dtype=float).reshape(3, 3) - np.eye(3)
+        for pair in usable
+        if pair.t_B_m is not None
+    ]
+    if translation_rows:
+        j_t = np.vstack(translation_rows)
+        singular_t = np.linalg.svd(j_t, compute_uv=False)
+        cond_t = float(singular_t[0] / max(singular_t[-1], 1e-12))
+        translation_information = float(
+            singular_t[-1] / np.sqrt(max(len(translation_rows), 1))
+        )
+    else:
+        cond_t = 1e9
+        translation_information = 0.0
+    translation_ok = (
+        len(translation_rows) >= 5
+        and cond_t < condition_threshold
+        and translation_information > 0.02
+    )
 
     if not rotation_ok:
-        notes.append(f"rotation condition {cond_r:.1f} exceeds threshold {condition_threshold}")
+        notes.append(
+            f"rotation not observable: condition={cond_r:.1f}, "
+            f"min_information={rotation_information:.3e}"
+        )
     else:
-        notes.append(f"rotation condition {cond_r:.1f}")
+        notes.append(
+            f"rotation observable: condition={cond_r:.1f}, "
+            f"min_information={rotation_information:.3e}"
+        )
     if not translation_ok:
         notes.append(
-            f"translation not observable (axis_rank={axis_rank}, trans_span={trans_span:.3f} m); "
-            "V1 will reject full SE3 without strong priors"
+            f"translation not observable: condition={cond_t:.1f}, "
+            f"min_information={translation_information:.3e}; "
+            "full SE3 will be rejected"
         )
     return ObservabilityReport(
         rotation_observable=rotation_ok,

imu_lidar/phase_a_replay.py

@@ -0,0 +1,290 @@
+"""Cached Phase-A replay: rehydrate gyro factors, compare variants, write reports."""
+
+from __future__ import annotations
+
+import json
+from collections import defaultdict
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+
+from .imu_io import load_imu_samples
+from .motion_pairs_io import (
+    build_motion_pairs_payload,
+    load_motion_pairs,
+    pair_from_dict,
+    save_motion_pairs,
+)
+from .phase_a import (
+    ProgressCallback,
+    phase_a_comparison_to_dict,
+    phase_a_metadata_complete,
+    rehydrate_phase_a_pairs,
+    solve_phase_a_comparison,
+)
+from .vehicle_config import load_vehicle_config, prior_enabled
+
+
+def _rotation_prior(
+    vehicle_config_path: Path,
+) -> tuple[np.ndarray | None, float]:
+    config = load_vehicle_config(vehicle_config_path)
+    if not prior_enabled(config, "rotation_prior"):
+        return None, 15.0
+    prior = (config.get("initialization") or {}).get("rotation_prior") or {}
+    matrix = prior.get("R_IMU_lidar")
+    if matrix is None:
+        return None, float(prior.get("sigma_deg", 15.0))
+    return (
+        np.asarray(matrix, dtype=float).reshape(3, 3),
+        float(prior.get("sigma_deg", 15.0)),
+    )
+
+
+def _sanitize_json(value: Any) -> Any:
+    if isinstance(value, dict):
+        return {str(key): _sanitize_json(item) for key, item in value.items()}
+    if isinstance(value, (list, tuple)):
+        return [_sanitize_json(item) for item in value]
+    if isinstance(value, np.ndarray):
+        return _sanitize_json(value.tolist())
+    if isinstance(value, (np.floating, float)):
+        number = float(value)
+        return number if np.isfinite(number) else None
+    if isinstance(value, (np.integer, np.bool_)):
+        return value.item()
+    return value
+
+
+def _write_json(path: Path, payload: Any) -> None:
+    path.write_text(
+        json.dumps(_sanitize_json(payload), indent=2, ensure_ascii=False) + "\n",
+        encoding="utf-8",
+    )
+
+
+def _load_cached_sessions(
+    motion_pairs_path: Path,
+) -> tuple[
+    dict[str, Any],
+    list,
+    dict[str, np.ndarray],
+    dict[str, float],
+]:
+    payload = load_motion_pairs(motion_pairs_path)
+    pairs = []
+    biases: dict[str, np.ndarray] = {}
+    offsets: dict[str, float] = {}
+    for session in payload.get("sessions") or []:
+        session_id = str(session["session_id"])
+        biases[session_id] = np.asarray(
+            session.get("gyro_bias_rad_s", np.zeros(3)),
+            dtype=float,
+        ).reshape(3)
+        offsets[session_id] = float(session.get("delta_t_s", 0.0))
+        pairs.extend(
+            pair_from_dict(item)
+            for item in session.get("pairs") or []
+        )
+    if not pairs:
+        raise ValueError(f"motion-pair cache is empty: {motion_pairs_path}")
+    return payload, pairs, biases, offsets
+
+
+def run_phase_a_replay(
+    *,
+    motion_pairs_path: Path,
+    vehicle_config_path: Path,
+    output_directory: Path,
+    imu_paths_by_session: dict[str, Path] | None = None,
+    excluded_sessions: set[str] | None = None,
+    strong_rotation_min_deg: float = 1.0,
+    decorrelation_block_s: float = 3.0,
+    max_pairs_per_block: int = 1,
+    bias_prior_sigma_rad_s: float = 0.002,
+    yaw_std_max_deg: float = 0.5,
+    leave_one_out_yaw_range_max_deg: float = 1.0,
+    data_prior_difference_max_deg: float = 1.0,
+    max_nfev: int = 200,
+    progress_callback: ProgressCallback | None = None,
+) -> dict[str, Any]:
+    """Run Phase-A only. Existing LiDAR relative motions are never recomputed."""
+
+    output_directory.mkdir(parents=True, exist_ok=True)
+    source_payload, pairs, bias0, offsets = _load_cached_sessions(
+        motion_pairs_path
+    )
+    session_ids = sorted(bias0)
+    if progress_callback is not None:
+        progress_callback(
+            "cache_loaded",
+            {
+                "schema_version": source_payload.get("schema_version"),
+                "sessions": len(session_ids),
+                "pairs": len(pairs),
+            },
+        )
+
+    rehydration_report: dict[str, Any] = {
+        "required": not phase_a_metadata_complete(pairs),
+        "pair_count": len(pairs),
+    }
+    if not phase_a_metadata_complete(pairs):
+        supplied_paths = {} if imu_paths_by_session is None else imu_paths_by_session
+        missing = [sid for sid in session_ids if sid not in supplied_paths]
+        if missing:
+            raise ValueError(
+                "v1 cache lacks J_bg/cov; provide --session-imu for: "
+                + ", ".join(missing)
+            )
+        imu_by_session = {
+            sid: load_imu_samples(supplied_paths[sid])
+            for sid in session_ids
+        }
+        pairs, details = rehydrate_phase_a_pairs(
+            pairs,
+            imu_by_session=imu_by_session,
+            bias0_by_session=bias0,
+            progress_callback=progress_callback,
+        )
+        rehydration_report.update(details)
+        if float(details["max_R_A_error_deg"]) > 0.05:
+            raise ValueError(
+                "rehydrated IMU rotations do not match cached R_A: "
+                f"max error={details['max_R_A_error_deg']:.6f} deg; "
+                "check session-to-IMU path mapping"
+            )
+
+    grouped: dict[str, list] = defaultdict(list)
+    for pair in pairs:
+        grouped[pair.session_id].append(pair)
+    enriched_payload = build_motion_pairs_payload(
+        prepared_sessions=[
+            {
+                "session_id": sid,
+                "time_offset_s": offsets[sid],
+                "gyro_bias_rad_s": bias0[sid],
+                "pairs": tuple(grouped[sid]),
+            }
+            for sid in session_ids
+        ]
+    )
+    enriched_cache_path = save_motion_pairs(
+        output_directory / "motion_pairs_phase_a_v2.json",
+        enriched_payload,
+    )
+
+    rotation_prior, rotation_prior_sigma_deg = _rotation_prior(
+        vehicle_config_path
+    )
+    comparison = solve_phase_a_comparison(
+        pairs,
+        gyro_bias_rad_s_by_session=bias0,
+        rotation_prior=rotation_prior,
+        rotation_prior_sigma_deg=rotation_prior_sigma_deg,
+        preexcluded_session_ids=excluded_sessions,
+        strong_rotation_min_deg=strong_rotation_min_deg,
+        decorrelation_block_s=decorrelation_block_s,
+        max_pairs_per_block=max_pairs_per_block,
+        bias_prior_sigma_rad_s=bias_prior_sigma_rad_s,
+        yaw_std_max_deg=yaw_std_max_deg,
+        leave_one_out_yaw_range_max_deg=(
+            leave_one_out_yaw_range_max_deg
+        ),
+        data_prior_difference_max_deg=data_prior_difference_max_deg,
+        run_leave_one_out=True,
+        max_nfev=max_nfev,
+        progress_callback=progress_callback,
+    )
+    full = phase_a_comparison_to_dict(comparison)
+    full["input"] = {
+        "motion_pairs": str(motion_pairs_path),
+        "source_schema_version": source_payload.get("schema_version"),
+        "vehicle_config": str(vehicle_config_path),
+        "session_imu_paths": {
+            sid: str(path)
+            for sid, path in (imu_paths_by_session or {}).items()
+        },
+        "excluded_sessions": sorted(excluded_sessions or set()),
+    }
+    full["rehydration"] = rehydration_report
+    full["enriched_cache"] = str(enriched_cache_path)
+    full["parameters"] = {
+        "strong_rotation_min_deg": strong_rotation_min_deg,
+        "decorrelation_block_s": decorrelation_block_s,
+        "max_pairs_per_block": max_pairs_per_block,
+        "bias_prior_sigma_rad_s": bias_prior_sigma_rad_s,
+        "rotation_prior_sigma_deg": rotation_prior_sigma_deg,
+        "yaw_std_max_deg": yaw_std_max_deg,
+        "leave_one_out_yaw_range_max_deg": (
+            leave_one_out_yaw_range_max_deg
+        ),
+        "data_prior_difference_max_deg": (
+            data_prior_difference_max_deg
+        ),
+        "max_nfev": max_nfev,
+    }
+
+    variants = full["variants"]
+    summary = {
+        "status": comparison.solution_status,
+        "accepted": comparison.accepted,
+        "partial_accepted": comparison.partial_accepted,
+        "acceptance_checks": comparison.acceptance_checks,
+        "primary_result": comparison.recommended_result,
+        "variants": {
+            name: {
+                "rpy_deg_xyz": item["rpy_deg_xyz"],
+                "R_IMU_lidar": item["R_IMU_lidar"],
+                "residual_rms_deg": item["residual_rms_deg"],
+                "residual_p95_deg": item["residual_p95_deg"],
+                "accepted": item["accepted"],
+                "gyro_bias_rad_s_per_session": item[
+                    "gyro_bias_rad_s_per_session"
+                ],
+            }
+            for name, item in variants.items()
+            if item is not None
+        },
+        "marginal_observability_A1": full[
+            "marginal_observability_A1"
+        ],
+        "data_vs_prior_yaw_diff_deg": (
+            comparison.data_vs_prior_yaw_diff_deg
+        ),
+        "data_vs_prior_geodesic_deg": (
+            comparison.data_vs_prior_geodesic_deg
+        ),
+        "leave_one_out_yaw_range_deg": (
+            comparison.leave_one_out_yaw_range_deg
+        ),
+        "leave_one_out_observable_max_deg": (
+            comparison.leave_one_out_observable_max_deg
+        ),
+        "strong_pair_candidate_count": (
+            comparison.strong_pair_candidate_count
+        ),
+        "decorrelated_pair_count": comparison.decorrelated_pair_count,
+        "strong_pair_counts_per_session": (
+            comparison.strong_pair_counts_per_session
+        ),
+        "excluded_sessions": list(comparison.excluded_sessions),
+        "rehydration": rehydration_report,
+        "comparison_file": "phase_a_comparison.json",
+        "observability_file": "phase_a_observability.json",
+        "leave_one_out_file": "phase_a_leave_one_out.json",
+        "enriched_cache_file": enriched_cache_path.name,
+    }
+
+    _write_json(output_directory / "phase_a_comparison.json", full)
+    _write_json(
+        output_directory / "phase_a_observability.json",
+        full["marginal_observability_A1"],
+    )
+    _write_json(
+        output_directory / "phase_a_leave_one_out.json",
+        full["leave_one_out"],
+    )
+    _write_json(output_directory / "phase_a_summary.json", summary)
+    return summary

imu_lidar/rotation_handeye.py

@@ -16,6 +16,8 @@ class RotationHandeyeResult:
     R_IMU_lidar: np.ndarray
     residual_rms_deg: float
     residual_median_deg: float
+    residual_p95_deg: float
+    outlier_fraction_gt_5deg: float
     pair_count: int
     ok: bool
     notes: tuple[str, ...] = ()
@@ -28,17 +30,21 @@ def _pair_weight(pair: MotionPair) -> float:
     return weight
 
 
-def _tsai_rotation_initial(pairs: list[MotionPair]) -> np.ndarray:
+def _tsai_rotation_initial(
+    pairs: list[MotionPair],
+    pair_weights: np.ndarray | None = None,
+) -> np.ndarray:
     """Closed-form rotation hand-eye initial guess (Tsai-style linear solve)."""
 
     rows: list[np.ndarray] = []
     rhs: list[np.ndarray] = []
-    for pair in pairs:
+    weights = np.ones(len(pairs)) if pair_weights is None else np.asarray(pair_weights, dtype=float)
+    for pair, pair_weight in zip(pairs, weights):
         alpha = so3_log(pair.R_A)
         beta = so3_log(pair.R_B)
         if np.linalg.norm(alpha) < 1e-6 or np.linalg.norm(beta) < 1e-6:
             continue
-        w = np.sqrt(_pair_weight(pair))
+        w = np.sqrt(float(pair_weight))
         rows.append(w * skew(alpha + beta))
         rhs.append(w * (beta - alpha))
     if len(rows) < 2:
@@ -57,26 +63,101 @@ def _pair_residual_deg(r_x: np.ndarray, pair: MotionPair) -> float:
     return float(np.degrees(np.linalg.norm(err)))
 
 
-def solve_rotation_handeye(pairs: list[MotionPair] | tuple[MotionPair, ...]) -> RotationHandeyeResult:
-    """Solve ``R_A R_X = R_X R_B`` with weighted robust nonlinear refinement."""
+def _rms_deg(r_x: np.ndarray, pairs: list[MotionPair]) -> float:
+    if not pairs:
+        return 1e9
+    errs = np.asarray([_pair_residual_deg(r_x, pair) for pair in pairs], dtype=float)
+    return float(np.sqrt(np.mean(errs**2)))
 
-    usable = [pair for pair in pairs if rotation_angle_deg(pair.R_A) > 1.0 and rotation_angle_deg(pair.R_B) > 1.0]
+
+def select_strong_rotation_pairs(
+    pairs: list[MotionPair] | tuple[MotionPair, ...],
+    *,
+    min_rotation_deg: float = 1.0,
+) -> list[MotionPair]:
+    """Return pairs that independently excite rotation on both sensor sides."""
+
+    threshold = float(min_rotation_deg)
+    return [
+        pair
+        for pair in pairs
+        if rotation_angle_deg(pair.R_A) > threshold
+        and rotation_angle_deg(pair.R_B) > threshold
+    ]
+
+
+def estimate_rotation_handeye_initial(
+    pairs: list[MotionPair] | tuple[MotionPair, ...],
+    *,
+    min_rotation_deg: float = 1.0,
+) -> np.ndarray:
+    """Return the fast data-only Tsai initialization without nonlinear refine."""
+
+    usable = select_strong_rotation_pairs(
+        pairs,
+        min_rotation_deg=min_rotation_deg,
+    )
+    if not usable:
+        return np.eye(3)
+    raw_weights = np.asarray(
+        [_pair_weight(pair) for pair in usable],
+        dtype=float,
+    )
+    median = max(float(np.median(raw_weights)), 1e-12)
+    weights = np.clip(raw_weights / median, 0.1, 10.0)
+    return _tsai_rotation_initial(usable, weights)
+
+
+def solve_rotation_handeye(
+    pairs: list[MotionPair] | tuple[MotionPair, ...],
+    *,
+    R_prior: np.ndarray | None = None,
+    prior_sigma_deg: float | None = None,
+) -> RotationHandeyeResult:
+    """Solve ``R_A R_X = R_X R_B`` with weighted robust nonlinear refinement.
+
+    Optional CAD / installation ``R_prior`` soft-constrains the extrinsic yaw that
+    is weakly observable under near-planar motion.
+    """
+
+    usable = select_strong_rotation_pairs(pairs)
     notes: list[str] = []
     if len(usable) < 3:
         return RotationHandeyeResult(
             R_IMU_lidar=np.eye(3),
             residual_rms_deg=1e9,
             residual_median_deg=1e9,
+            residual_p95_deg=1e9,
+            outlier_fraction_gt_5deg=1.0,
             pair_count=len(usable),
             ok=False,
             notes=("need at least 3 motion pairs with meaningful rotation",),
         )
 
-    r0 = _tsai_rotation_initial(usable)
-    weights = np.asarray([_pair_weight(pair) for pair in usable], dtype=float)
+    raw_weights = np.asarray([_pair_weight(pair) for pair in usable], dtype=float)
+    median_raw_weight = max(float(np.median(raw_weights)), 1e-12)
+    weights = np.clip(raw_weights / median_raw_weight, 0.1, 10.0)
+    r0 = _tsai_rotation_initial(usable, weights)
+    r_prior = None
+    if R_prior is not None:
+        r_prior = orthonormalize_rotation(np.asarray(R_prior, dtype=float).reshape(3, 3))
+        rms_tsai = _rms_deg(r0, usable)
+        rms_prior = _rms_deg(r_prior, usable)
+        if rms_prior <= rms_tsai * 1.25:
+            r0 = r_prior
+            notes.append(
+                f"init from rotation prior (rms={rms_prior:.3f} deg vs Tsai {rms_tsai:.3f} deg)"
+            )
+        else:
+            notes.append(
+                f"init from Tsai (rms={rms_tsai:.3f} deg; prior {rms_prior:.3f} deg kept as soft constraint)"
+            )
+
     notes.append(
-        f"weighted hand-eye: weight median={float(np.median(weights)):.3g}, "
-        f"min={float(np.min(weights)):.3g}, max={float(np.max(weights)):.3g}"
+        "weighted hand-eye: normalized/clipped IMU confidence "
+        f"raw_median={median_raw_weight:.3g}, "
+        f"normalized_min={float(np.min(weights)):.3g}, "
+        f"normalized_max={float(np.max(weights)):.3g}"
     )
 
     def pack(r: np.ndarray) -> np.ndarray:
@@ -85,12 +166,21 @@ def solve_rotation_handeye(pairs: list[MotionPair] | tuple[MotionPair, ...]) ->
     def unpack(vec: np.ndarray) -> np.ndarray:
         return orthonormalize_rotation(so3_exp(vec))
 
+    sigma = 15.0 if prior_sigma_deg is None else float(prior_sigma_deg)
+    prior_w = 0.0
+    if r_prior is not None and sigma > 1e-6:
+        # Scale prior to a few strong pairs so it regularizes yaw without dominating.
+        prior_w = float(np.sqrt(np.median(weights)) / np.deg2rad(sigma))
+        notes.append(f"rotation prior soft constraint sigma={sigma:.1f} deg, weight={prior_w:.3g}")
+
     def residual(vec: np.ndarray) -> np.ndarray:
         r_x = unpack(vec)
         residuals = []
         for pair, weight in zip(usable, weights):
             err = so3_log(r_x.T @ pair.R_A @ r_x @ pair.R_B.T)
             residuals.append(np.sqrt(weight) * err)
+        if r_prior is not None and prior_w > 0:
+            residuals.append(prior_w * so3_log(r_prior.T @ r_x))
         return np.concatenate(residuals)
 
     opt = least_squares(residual, pack(r0), loss="huber", f_scale=np.deg2rad(1.0), max_nfev=200)
@@ -99,14 +189,28 @@ def solve_rotation_handeye(pairs: list[MotionPair] | tuple[MotionPair, ...]) ->
     # Report unweighted RMS/median for interpretability.
     rms = float(np.sqrt(np.mean(errs**2)))
     med = float(np.median(errs))
+    p95 = float(np.percentile(errs, 95.0))
+    outlier_fraction = float(np.mean(errs > 5.0))
     notes.append(f"optimized over {len(usable)} pairs")
-    ok = rms < 5.0 and len(usable) >= 3
+    notes.append(
+        f"rotation residual quality: rms={rms:.3f} deg, median={med:.3f} deg, "
+        f"p95={p95:.3f} deg, >5deg={100.0 * outlier_fraction:.2f}%"
+    )
+    ok = (
+        len(usable) >= 3
+        and rms < 1.5
+        and med < 0.5
+        and p95 < 1.5
+        and outlier_fraction <= 0.005
+    )
     if not ok:
-        notes.append("rotation residual RMS too high or too few pairs")
+        notes.append("rotation residual distribution failed acceptance gates")
     return RotationHandeyeResult(
         R_IMU_lidar=r_x,
         residual_rms_deg=rms,
         residual_median_deg=med,
+        residual_p95_deg=p95,
+        outlier_fraction_gt_5deg=outlier_fraction,
         pair_count=len(usable),
         ok=ok,
         notes=tuple(notes),

imu_lidar/time_offset.py

@@ -73,7 +73,10 @@ def _correlate_offset(
     y0, y1, y2 = peaks
     denom = y0 - 2 * y1 + y2
     if abs(denom) > 1e-12:
-        best_delta = float(best_delta + 0.5 * (y0 - y2) / denom * dt)
+        refined = float(best_delta + 0.5 * (y0 - y2) / denom * dt)
+        # Parabola can jump outside the searched window; keep it clamped.
+        if abs(refined) <= search_s + dt:
+            best_delta = refined
         best_peak = float(y1)
     return best_delta, best_peak
 
@@ -99,9 +102,12 @@ def estimate_time_offset(
     bias = np.zeros(3) if gyro_bias_rad_s is None else np.asarray(gyro_bias_rad_s, dtype=float)
     gyro = imu.gyro_rad_s - bias
 
-    stride = max(1, len(frames) // 20)
+    # Use short consecutive (or near-consecutive) pairs. A large stride (e.g.
+    # len//20) averages over many seconds and destroys |ω| correlation even when
+    # host/device clocks are already aligned.
+    stride = 1 if len(frames) < 80 else 2
     rotations, pair_times = estimate_frame_rotations(frames, stride=stride)
-    if len(rotations) < 4:
+    if len(rotations) < 8:
         rotations, pair_times = estimate_frame_rotations(frames, stride=1)
     if len(rotations) < 4:
         return TimeOffsetResult(0.0, 0.0, search_s, ("not enough LiDAR relative rotations",), False)
@@ -131,8 +137,16 @@ def estimate_time_offset(
         f"LiDAR mean pair rotation {np.mean([rotation_angle_deg(r) for r in rotations]):.2f} deg"
     )
     notes.append(f"searched delta_t in ±{search_s:.3f}s by direct correlation")
-    ok = peak > 0.15
-    if not ok:
+    # Host-UTC-bridged sessions are already on one timeline; |ω| peak can stay
+    # weak even at the correct lag (ICP rate vs gyro scale). Accept near-zero δt.
+    near_zero = abs(float(delta)) <= min(0.05, 0.25 * float(search_s))
+    ok = peak > 0.15 or near_zero
+    if peak <= 0.15 and near_zero:
+        notes.append(
+            f"correlation peak weak ({peak:.3f}) but |delta_t|={abs(delta):.4f}s ~0; "
+            "accepting as already-aligned (e.g. host UTC bridge)"
+        )
+    elif not ok:
         notes.append("correlation peak is weak; check overlapping motion and axis units")
     return TimeOffsetResult(
         delta_t_s=delta,
@@ -176,6 +190,7 @@ def refine_time_offset_signed(
     gyro_bias_rad_s: np.ndarray | None = None,
     search_s: float = 0.08,
     sample_hz: float = 50.0,
+    max_shift_s: float | None = 0.05,
 ) -> TimeOffsetResult:
     """Refine ``δt`` with signed 3-axis rates using a known ``R_IMU_lidar``.
 
@@ -279,9 +294,23 @@ def refine_time_offset_signed(
         f"corr={best_corr:.3f}, mag_corr={mag_at_best:.3f} (coarse_mag={mag_at_coarse:.3f}), "
         f"search=±{half:.3f}s"
     )
+    shift = abs(best_delta - float(delta_t_s))
+    if max_shift_s is not None and shift > float(max_shift_s):
+        notes.append(
+            f"signed refine rejected: |Δδt|={shift:.4f}s exceeds max_shift={float(max_shift_s):.4f}s; "
+            "keeping previous delta_t"
+        )
+        return TimeOffsetResult(
+            delta_t_s=float(delta_t_s),
+            correlation_peak=mag_at_coarse if mag_at_coarse > 0 else best_corr,
+            search_s=search_s,
+            notes=tuple(notes),
+            ok=True,
+        )
+    # Require a meaningful MSE drop so tiny downhill noise cannot walk δt across iterations.
     improved = (
         np.isfinite(best_cost)
-        and best_cost < coarse_cost * 0.999
+        and best_cost < coarse_cost * 0.98
         # Do not sacrifice the more reliable magnitude alignment for a noisy signed MSE gain.
         and mag_at_best + 1e-4 >= mag_at_coarse
     )

imu_lidar/文件职责说明.md

@@ -41,10 +41,11 @@ python -m imu_lidar.cli run --vehicle-config ... --imu ... --lidar ... --output
 | 5   | `lidar_deskew.py`       | 可选点云去畸变（低速可关）                 |
 | 6   | `imu_preintegration.py` | IMU 预积分（旋转及速度/位移增量、协方差、零偏雅可比） |
 | 6   | `motion_pairs.py`       | 构造运动对；手眼使用其中的旋转               |
+| 6   | `motion_pairs_io.py`    | 运动对 JSON 缓存读写（供可视化直读）          |
 | 7   | `rotation_handeye.py`   | 加权旋转手眼                        |
 | 8   | `observability.py`      | 旋转 / 平移可观性检查                  |
 | 8   | `joint_optimizer.py`    | 联合精修；完整模式下可估计平移、重力、速度与时变零偏    |
-| 9   | `finalize.py`           | 写出结果 JSON                     |
+| 9   | `finalize.py`           | 写出结果 JSON（含 `motion_pairs.json`） |
 | —   | `pipeline.py`           | 编排全流程                         |
 | —   | `cli.py`                | 命令行入口                         |
 | —   | `CHANGELOG.md`          | 改动记录                          |
@@ -56,7 +57,7 @@ python -m imu_lidar.cli run --vehicle-config ... --imu ... --lidar ... --output
 
 - **运动对**始终计算完整预积分量（旋转、速度增量、位移增量及不确定度）。
 - `--mode rotation_only`：只精修旋转与常值陀螺零偏，交付旋转与时间偏置。
-- `--mode full_se3`：在可观时再估计重力、关键帧速度、时变零偏与平移；结果写入 `summary.json` 的 joint 字段。
+- `--mode full_se3`：当前完成 Phase-A 后明确拒绝平移；待 Phase-B/C 会话状态重构完成后再恢复完整 SE(3) 交付。
 
 ---
 

tests/README.md

@@ -28,8 +28,14 @@ python -m pytest -q
 | `test_preintegration_bias_jacobian_matches_finite_difference` | 随机陀螺序列       | 旋转预积分 `J_bg`        | 与有限差分一致（松阈值）                               |
 | `test_imu_preintegration_recovers_constant_accel_translation` | 常值加速度        | 完整预积分 Δv/Δp         | 接近解析值                                      |
 | `test_imu_preintegration_bias_jacobian_finite_difference`     | 随机 IMU       | `J_bg`/`J_ba` 一阶修正  | 与重积分接近                                     |
-| `test_synthetic_pipeline_rotation_and_time_offset`            | 端到端合成会话      | `rotation_only` 全流程 | `rotation_only_accepted`；δt 准；手眼 RMS < 5°  |
-| `test_synthetic_pipeline_full_se3_smoke`                      | 同上           | `full_se3` 不崩溃      | 状态为 accepted / rejected / rotation_only 之一 |
+| `test_synthetic_pipeline_rejects_noisy_icp_but_keeps_time_audit` | synthetic end-to-end | strict rotation quality gate + time audit | noisy ICP is blocked; delta-t remains accurate |
+| `test_synthetic_pipeline_full_se3_smoke` | synthetic end-to-end | full-SE(3) smoke test | returns an explicit accepted/rejected/blocked status |
+| `test_planar_yaw_is_not_full_rotation_or_translation_observable` | pure-yaw motion pairs | degeneracy detection | full rotation/translation observability is rejected |
+| `test_multi_axis_motion_is_rotation_and_translation_observable` | multi-axis motion pairs | positive observability case | rotation and translation pass |
+| `test_translation_prior_is_reported_but_not_accepted_when_unobservable` | planar motion + CAD prior | prior semantics | prior is reported but not accepted as calibration |
+| `test_handeye_rejects_a_small_fraction_of_gross_rotation_outliers` | motion pairs with a gross outlier | residual-distribution gate | solve is rejected |
+| `test_motion_pairs_reject_low_fitness` | low-fitness registration | fitness gate | no motion pair is emitted |
+| `test_motion_pairs_reject_imu_and_lidar_discontinuities` | timestamp gaps | continuity gates | cross-gap pairs are rejected |
 
 
 合成数据由 `tools/generate_synthetic_session.py` 生成（墙面点云 + 已知外参 yaw 与 δt）。  

tests/test_h32_dlog_export.py

@@ -8,7 +8,8 @@ from pathlib import Path
 import numpy as np
 
 from tools.h32_dlog.difop import CHANNELS, HORIZONTAL_START, VERTICAL_START, parse_difop_angles
-from tools.h32_dlog.dobject import discover_records, iter_payloads, resolve_dlog_root
+from tools.h32_dlog.dobject import RECORD_RE, discover_records, iter_payloads, resolve_dlog_root
+from tools.h32_dlog.timeutil import local_wall_to_dotnet_ticks
 from tools.h32_dlog.load_session import load_h32_dlog_lidar
 from tools.h32_dlog.payload_v1 import (
     MsopPacketItem,
@@ -90,6 +91,19 @@ def _write_dorec_record(
     return start
 
 
+def test_recovered_index_line_and_local_ticks():
+    line = (
+        ">DObject `frontlidar-msop-raw` post len=15532B, id:9CF1, "
+        "tic:639218060782100466, @data.bin:0"
+    )
+    match = RECORD_RE.search(line)
+    assert match is not None
+    assert match.group("name") == "frontlidar-msop-raw"
+    assert match.group("file") == "data.bin"
+    assert int(match.group("offset")) == 0
+    assert local_wall_to_dotnet_ticks("2026-08-08T17:14:38") == 639218060780000000
+
+
 def test_parse_msop_and_difop_payload_roundtrip():
     packet = _make_msop_packet(seconds=1700000000, microseconds=123456)
     item = MsopPacketItem(

tests/test_hi13_imu.py

@@ -0,0 +1,79 @@
+"""Unit tests for HI13 / HI91 IMU decoding."""
+
+from __future__ import annotations
+
+import struct
+
+from tools.rscap_v2.capture_format_v2 import CaptureFile, CaptureHeader, RawChunk
+from tools.rscap_v2.hi13_imu import crc16_hi13, iter_hi13_imu_samples, parse_hi91_frame
+
+
+def _hi91_frame(
+    *,
+    device_ms: int = 123456,
+    accel_g=(0.0, 0.0, 1.0),
+    gyro_dps=(1.0, -2.0, 3.0),
+) -> bytes:
+    payload = bytearray(76)
+    payload[0] = 0x91
+    struct.pack_into("<H", payload, 1, 0)  # pps
+    payload[3] = 25  # temp
+    struct.pack_into("<f", payload, 4, 101325.0)
+    struct.pack_into("<I", payload, 8, device_ms)
+    struct.pack_into("<fff", payload, 12, *accel_g)
+    struct.pack_into("<fff", payload, 24, *gyro_dps)
+    # remaining mag/rpy/quat left zero
+    payload_length = len(payload)
+    header = bytearray(6)
+    header[0] = 0x5A
+    header[1] = 0xA5
+    header[2] = payload_length & 0xFF
+    header[3] = (payload_length >> 8) & 0xFF
+    frame_wo_crc = bytes(header[:4]) + bytes(payload)
+    # crc over header[0:4] + payload
+    tmp = bytearray(6 + payload_length)
+    tmp[0:4] = header[0:4]
+    tmp[6:] = payload
+    crc = crc16_hi13(tmp, payload_length)
+    header[4] = crc & 0xFF
+    header[5] = (crc >> 8) & 0xFF
+    return bytes(header) + bytes(payload)
+
+
+def test_parse_hi91_units():
+    frame = _hi91_frame(device_ms=5000, accel_g=(0.0, 0.0, 1.0), gyro_dps=(57.2957795, 0.0, 0.0))
+    parsed = parse_hi91_frame(frame)
+    assert parsed is not None
+    gyro, accel, device_ms = parsed
+    assert device_ms == 5000
+    assert abs(accel[2] - 9.80665) < 1e-4
+    assert abs(gyro[0] - 1.0) < 1e-5
+
+
+def test_iter_hi13_from_capture():
+    frame = _hi91_frame(device_ms=42)
+    header = CaptureHeader(
+        sensor_kind="hi13r4-imu",
+        session_id="t",
+        session_start_utc_ticks=0,
+        session_start_monotonic_ticks=0,
+        monotonic_frequency=10_000_000,
+        port="COM1",
+        baud=115200,
+        file_start_utc_ticks=0,
+    )
+    chunk = RawChunk(
+        sequence=1,
+        receive_utc_ticks=100,
+        receive_monotonic_ticks=1,
+        raw=frame,
+        record_file_offset=0,
+        raw_file_offset=0,
+        record_crc32=0,
+        crc_valid=True,
+    )
+    capture = CaptureFile(path="mem", header=header, chunks=[chunk], footer=None)
+    samples = iter_hi13_imu_samples(capture)
+    assert len(samples) == 1
+    assert samples[0].device_timestamp_us == 42_000
+    assert abs(samples[0].t_s - 0.042) < 1e-12

tests/test_motion_pairs_io.py

@@ -0,0 +1,78 @@
+"""Tests for motion-pair cache IO."""
+
+from __future__ import annotations
+
+from pathlib import Path
+
+import numpy as np
+
+from imu_lidar.contracts import MotionPair
+from imu_lidar.motion_pairs_io import (
+    build_motion_pairs_payload,
+    load_motion_pairs,
+    pair_from_dict,
+    pair_to_dict,
+    pairs_for_session,
+    save_motion_pairs,
+)
+
+
+def test_pair_roundtrip(tmp_path: Path) -> None:
+    pair = MotionPair(
+        session_id="s0",
+        i=1,
+        j=4,
+        t_i_s=1.0,
+        t_j_s=2.5,
+        R_A=np.eye(3),
+        R_B=np.eye(3),
+        t_A_m=np.array([0.1, 0.0, 0.0]),
+        t_B_m=np.array([0.1, 0.0, 0.0]),
+        fitness=0.8,
+        metadata={
+            "weight": 12.0,
+            "cov": (np.eye(3) * 1e-4).tolist(),
+            "J_bg": (-np.eye(3)).tolist(),
+            "cov9": [[0.0] * 9] * 9,
+            "backend": "test",
+            "gyro_bias0_rad_s": [0.01, -0.02, 0.03],
+            "accel_bias0_m_s2": [0.1, 0.2, -0.1],
+            "time_offset_s": 0.004,
+            "keyframe_span": 3,
+            "is_consecutive": False,
+        },
+    )
+    encoded = pair_to_dict(pair)
+    assert "cov9" not in encoded["metadata"]
+    assert "cov" in encoded["metadata"]
+    assert "J_bg" in encoded["metadata"]
+    assert encoded["metadata"]["weight"] == 12.0
+    restored = pair_from_dict(encoded)
+    assert restored.i == 1 and restored.j == 4
+    np.testing.assert_allclose(restored.t_A_m, [0.1, 0.0, 0.0])
+    np.testing.assert_allclose(restored.metadata["gyro_bias0_rad_s"], [0.01, -0.02, 0.03])
+    assert restored.metadata["keyframe_span"] == 3
+
+    payload = build_motion_pairs_payload(
+        prepared_sessions=[
+            {
+                "session_id": "s0",
+                "time_offset_s": 0.0,
+                "gyro_bias_rad_s": np.zeros(3),
+                "pairs": (pair,),
+            }
+        ]
+    )
+    path = save_motion_pairs(tmp_path / "motion_pairs.json", payload)
+    loaded = load_motion_pairs(path)
+    assert loaded["schema_version"] == 2
+    pairs = pairs_for_session(loaded, "s0")
+    assert len(pairs) == 1
+    assert pairs[0].session_id == "s0"
+
+    payload["schema_version"] = 1
+    legacy_path = save_motion_pairs(
+        tmp_path / "motion_pairs_v1.json", payload
+    )
+    legacy = load_motion_pairs(legacy_path)
+    assert legacy["schema_version"] == 1

tests/test_phase_a.py

@@ -0,0 +1,194 @@
+"""Tests for cached, session-balanced Phase-A comparison."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from imu_lidar.contracts import ImuSeries, MotionPair
+from imu_lidar.geometry import so3_exp, so3_log
+from imu_lidar.imu_preintegration import preintegrate_gyro
+from imu_lidar.phase_a import (
+    rehydrate_phase_a_pairs,
+    select_decorrelated_phase_a_pairs,
+    solve_phase_a_comparison,
+)
+
+
+def _phase_a_pair(
+    session_id: str,
+    index: int,
+    r_true: np.ndarray,
+    vector_deg: tuple[float, float, float],
+    bias0: np.ndarray,
+) -> MotionPair:
+    r_b = so3_exp(np.deg2rad(np.asarray(vector_deg, dtype=float)))
+    return MotionPair(
+        session_id=session_id,
+        i=index,
+        j=index + 1,
+        t_i_s=float(index),
+        t_j_s=float(index + 1),
+        R_A=r_true @ r_b @ r_true.T,
+        R_B=r_b,
+        t_A_m=np.zeros(3),
+        t_B_m=np.zeros(3),
+        metadata={
+            "J_bg": (-np.eye(3)).tolist(),
+            "cov": (np.eye(3) * 1e-5).tolist(),
+            "gyro_bias0_rad_s": bias0.tolist(),
+        },
+    )
+
+
+def test_phase_a_reports_three_variants_and_leave_one_session_out() -> None:
+    r_true = so3_exp(np.deg2rad(np.array([3.0, -2.0, 25.0])))
+    prior = so3_exp(np.deg2rad(np.array([0.0, 0.0, 0.2]))) @ r_true
+    vectors = (
+        (12.0, 0.0, 0.0),
+        (0.0, 15.0, 0.0),
+        (0.0, 0.0, 18.0),
+        (10.0, 8.0, 0.0),
+        (0.0, 11.0, 9.0),
+        (7.0, 0.0, 13.0),
+        (9.0, -5.0, 6.0),
+        (-6.0, 8.0, 11.0),
+        (5.0, 7.0, -9.0),
+    )
+    biases = {
+        "s0": np.array([0.001, -0.0005, 0.0002]),
+        "s1": np.array([-0.0004, 0.0008, -0.0001]),
+        "s2": np.array([0.0002, 0.0001, -0.0006]),
+    }
+    pairs: list[MotionPair] = []
+    index = 0
+    for sid, count in (("s0", 18), ("s1", 9), ("s2", 6)):
+        for local_index in range(count):
+            pairs.append(
+                _phase_a_pair(
+                    sid,
+                    index,
+                    r_true,
+                    vectors[local_index % len(vectors)],
+                    biases[sid],
+                )
+            )
+            index += 1
+
+    result = solve_phase_a_comparison(
+        pairs,
+        gyro_bias_rad_s_by_session=biases,
+        rotation_prior=prior,
+        rotation_prior_sigma_deg=15.0,
+        yaw_std_max_deg=1.0,
+        leave_one_out_yaw_range_max_deg=1.0,
+        data_prior_difference_max_deg=1.0,
+        decorrelation_block_s=0.0,
+        max_nfev=80,
+    )
+
+    assert result.accepted
+    assert result.strong_pair_counts_per_session == {
+        "s0": 18,
+        "s1": 9,
+        "s2": 6,
+    }
+    assert len(result.leave_one_out) == 3
+    assert result.marginal_observability.rank == 3
+    assert result.leave_one_out_yaw_range_deg < 0.1
+    for variant in (
+        result.fixed_bg_data_only,
+        result.session_bg_data_only,
+        result.session_bg_with_rotation_prior,
+    ):
+        error_deg = np.degrees(
+            np.linalg.norm(
+                so3_log(r_true.T @ variant.R_IMU_lidar)
+            )
+        )
+        assert error_deg < 0.1
+
+
+def test_rehydrate_phase_a_pairs_recovers_jacobian_without_lidar() -> None:
+    t_s = np.linspace(0.0, 1.0, 201)
+    gyro = np.tile(np.array([0.12, -0.04, 0.2]), (t_s.size, 1))
+    bias0 = np.array([0.01, -0.005, 0.002])
+    imu = ImuSeries(
+        t_s=t_s,
+        gyro_rad_s=gyro,
+        acc_m_s2=np.zeros((t_s.size, 3)),
+    )
+    preint = preintegrate_gyro(t_s, gyro, 0.1, 0.8, bias0)
+    pair = MotionPair(
+        session_id="s0",
+        i=0,
+        j=1,
+        t_i_s=0.1,
+        t_j_s=0.8,
+        R_A=preint.delta_R,
+        R_B=preint.delta_R,
+        metadata={
+            "t_i_imu_s": 0.1,
+            "t_j_imu_s": 0.8,
+            "gyro_bias0_rad_s": bias0.tolist(),
+            "preint_sigma_rad": preint.sigma_rad,
+        },
+    )
+
+    enriched, report = rehydrate_phase_a_pairs(
+        [pair],
+        imu_by_session={"s0": imu},
+        bias0_by_session={"s0": bias0},
+    )
+
+    assert "J_bg" in enriched[0].metadata
+    assert "cov" in enriched[0].metadata
+    assert report["max_R_A_error_deg"] < 1e-8
+
+
+def test_phase_a_time_blocks_do_not_count_overlapping_pairs_as_independent() -> None:
+    r_true = so3_exp(np.deg2rad(np.array([1.0, -2.0, 20.0])))
+    bias = np.zeros(3)
+    pairs = [
+        _phase_a_pair("s0", index, r_true, (5.0 + index, 2.0, 1.0), bias)
+        for index in range(9)
+    ]
+    selected = select_decorrelated_phase_a_pairs(
+        pairs,
+        block_s=3.0,
+        max_pairs_per_block=1,
+    )
+    assert len(selected) == 3
+    assert all(pair in pairs for pair in selected)
+
+
+def test_phase_a_planar_motion_is_partial_and_keeps_weak_direction_from_prior() -> None:
+    r_true = so3_exp(np.deg2rad(np.array([4.0, -3.0, 31.0])))
+    prior = so3_exp(np.deg2rad(np.array([0.2, -0.1, 0.4]))) @ r_true
+    biases = {"s0": np.zeros(3), "s1": np.zeros(3)}
+    pairs: list[MotionPair] = []
+    for session_index, sid in enumerate(biases):
+        for index in range(12):
+            pairs.append(
+                _phase_a_pair(
+                    sid,
+                    session_index * 100 + index,
+                    r_true,
+                    (0.0, 0.0, 8.0 + index),
+                    biases[sid],
+                )
+            )
+    result = solve_phase_a_comparison(
+        pairs,
+        gyro_bias_rad_s_by_session=biases,
+        rotation_prior=prior,
+        decorrelation_block_s=0.0,
+        yaw_std_max_deg=0.5,
+        run_leave_one_out=False,
+        max_nfev=80,
+    )
+    assert not result.accepted
+    assert result.partial_accepted
+    assert result.solution_status == "phase_a_partial_accepted"
+    assert result.marginal_observability.precision_rank == 2
+    assert result.observable_subspace_with_prior is not None
+    assert np.isinf(result.marginal_observability.direction_std_deg[0])

tests/test_quality_gates.py

@@ -0,0 +1,260 @@
+"""Regression tests for calibration quality, continuity, and observability gates."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from imu_lidar.contracts import ImuSeries, LidarFrame, MotionPair
+from imu_lidar.geometry import make_transform, so3_exp, so3_log
+from imu_lidar.joint_optimizer import solve_joint_extrinsic
+from imu_lidar.motion_pairs import build_motion_pairs
+from imu_lidar.observability import analyze_observability
+from imu_lidar.registration import RegistrationResult
+from imu_lidar.rotation_handeye import solve_rotation_handeye
+
+
+def _motion_pair(index: int, rotation_vector: np.ndarray) -> MotionPair:
+    rotation = so3_exp(np.asarray(rotation_vector, dtype=float))
+    return MotionPair(
+        session_id="synthetic",
+        i=index,
+        j=index + 1,
+        t_i_s=float(index),
+        t_j_s=float(index + 1),
+        R_A=rotation,
+        R_B=rotation,
+        t_A_m=np.zeros(3),
+        t_B_m=np.array([0.1, -0.03, 0.0]),
+        fitness=0.9,
+        metadata={
+            "J_bg": (-np.eye(3)).tolist(),
+            "cov": (np.eye(3) * 1e-4).tolist(),
+            "gyro_bias0_rad_s": [0.0, 0.0, 0.0],
+        },
+    )
+
+
+def _frame(frame_id: str, mid_s: float) -> LidarFrame:
+    return LidarFrame(
+        frame_id=frame_id,
+        t_start_s=mid_s - 0.01,
+        t_end_s=mid_s + 0.01,
+        points_xyz=np.zeros((64, 3)),
+    )
+
+
+def _registration(*, fitness: float = 0.9) -> RegistrationResult:
+    rotation = so3_exp(np.deg2rad(np.array([0.0, 0.0, 10.0])))
+    return RegistrationResult(
+        transform=make_transform(np.array([0.4, 0.0, 0.0]), rotation),
+        fitness=fitness,
+        rotation_deg=10.0,
+        translation_m=0.4,
+        backend="test",
+        ok=True,
+    )
+
+
+def test_planar_yaw_is_not_full_rotation_or_translation_observable():
+    pairs = [
+        _motion_pair(i, np.deg2rad(np.array([0.0, 0.0, angle_deg])))
+        for i, angle_deg in enumerate((5.0, 8.0, 12.0, 17.0, 23.0, 31.0))
+    ]
+
+    report = analyze_observability(pairs, np.eye(3))
+
+    assert not report.rotation_observable
+    assert not report.translation_observable
+
+
+def test_multi_axis_motion_is_rotation_and_translation_observable():
+    vectors_deg = (
+        (12.0, 0.0, 0.0),
+        (0.0, 15.0, 0.0),
+        (0.0, 0.0, 18.0),
+        (10.0, 8.0, 0.0),
+        (0.0, 11.0, 9.0),
+        (7.0, 0.0, 13.0),
+    )
+    pairs = [
+        _motion_pair(i, np.deg2rad(np.asarray(vector_deg)))
+        for i, vector_deg in enumerate(vectors_deg)
+    ]
+
+    report = analyze_observability(pairs, np.eye(3))
+
+    assert report.rotation_observable
+    assert report.translation_observable
+
+
+def test_translation_prior_is_reported_but_not_accepted_when_unobservable():
+    pairs = [
+        _motion_pair(i, np.deg2rad(np.array([0.0, 0.0, angle_deg])))
+        for i, angle_deg in enumerate((5.0, 8.0, 12.0, 17.0, 23.0, 31.0))
+    ]
+    prior = np.array([0.3, -0.2, 0.5])
+
+    result = solve_joint_extrinsic(
+        pairs,
+        np.eye(3),
+        force_rotation_only=False,
+        enable_phase_c=False,
+        t_prior_m=prior,
+    )
+
+    assert not result.translation_accepted
+    np.testing.assert_allclose(result.T_IMU_lidar[:3, 3], prior)
+    assert any("prior only" in note for note in result.notes)
+
+
+def test_handeye_rejects_a_small_fraction_of_gross_rotation_outliers():
+    rng = np.random.default_rng(7)
+    r_true = so3_exp(np.deg2rad(np.array([2.0, -3.0, 20.0])))
+    pairs: list[MotionPair] = []
+    for index in range(100):
+        axis = rng.normal(size=3)
+        axis /= np.linalg.norm(axis)
+        r_b = so3_exp(axis * np.deg2rad(rng.uniform(8.0, 30.0)))
+        r_a = r_true @ r_b @ r_true.T
+        if index == 0:
+            r_a = so3_exp(np.deg2rad(np.array([18.0, 0.0, 0.0]))) @ r_a
+        pairs.append(
+            MotionPair(
+                session_id="outlier",
+                i=index,
+                j=index + 1,
+                t_i_s=float(index),
+                t_j_s=float(index + 1),
+                R_A=r_a,
+                R_B=r_b,
+            )
+        )
+
+    result = solve_rotation_handeye(pairs)
+
+    assert not result.ok
+    assert result.outlier_fraction_gt_5deg > 0.005
+
+
+def test_motion_pairs_reject_low_fitness(monkeypatch):
+    monkeypatch.setattr(
+        "imu_lidar.motion_pairs.register_lidar_pair",
+        lambda *_args, **_kwargs: _registration(fitness=0.3),
+    )
+    imu = ImuSeries(
+        t_s=np.linspace(0.0, 1.2, 121),
+        gyro_rad_s=np.zeros((121, 3)),
+        acc_m_s2=np.zeros((121, 3)),
+    )
+
+    result = build_motion_pairs(
+        session_id="fitness",
+        keyframes=[_frame("0", 0.1), _frame("1", 1.1)],
+        keyframe_indices=[0, 1],
+        imu=imu,
+        delta_t_s=0.0,
+        min_registration_fitness=0.5,
+    )
+
+    assert not result.pairs
+    assert any("fitness<0.50: 1" in note for note in result.notes)
+
+
+def test_motion_pairs_reject_imu_and_lidar_discontinuities(monkeypatch):
+    monkeypatch.setattr(
+        "imu_lidar.motion_pairs.register_lidar_pair",
+        lambda *_args, **_kwargs: _registration(),
+    )
+    imu_with_gap = ImuSeries(
+        t_s=np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.8, 0.9, 1.0, 1.1, 1.2]),
+        gyro_rad_s=np.zeros((10, 3)),
+        acc_m_s2=np.zeros((10, 3)),
+    )
+    imu_result = build_motion_pairs(
+        session_id="imu-gap",
+        keyframes=[_frame("0", 0.1), _frame("1", 1.1)],
+        keyframe_indices=[0, 1],
+        imu=imu_with_gap,
+        delta_t_s=0.0,
+        max_imu_gap_s=0.2,
+    )
+
+    assert not imu_result.pairs
+    assert any("IMU gap>0.200s: 1" in note for note in imu_result.notes)
+
+    continuous_imu = ImuSeries(
+        t_s=np.linspace(0.0, 2.2, 221),
+        gyro_rad_s=np.zeros((221, 3)),
+        acc_m_s2=np.zeros((221, 3)),
+    )
+    lidar_result = build_motion_pairs(
+        session_id="lidar-gap",
+        keyframes=[_frame("0", 0.1), _frame("2", 2.1)],
+        keyframe_indices=[0, 2],
+        imu=continuous_imu,
+        delta_t_s=0.0,
+        all_frame_times_s=np.array([0.1, 0.2, 2.1]),
+        max_lidar_gap_s=0.5,
+    )
+
+    assert not lidar_result.pairs
+    assert any("LiDAR gap>0.500s: 1" in note for note in lidar_result.notes)
+
+def test_phase_a_keeps_session_bias_linearization_points_independent():
+    r_true = so3_exp(np.deg2rad(np.array([2.0, -3.0, 20.0])))
+    bias0_by_session = {
+        "s0": np.array([0.010, -0.004, 0.002]),
+        "s1": np.array([-0.006, 0.008, -0.003]),
+    }
+    vectors_deg = (
+        (12.0, 0.0, 0.0),
+        (0.0, 15.0, 0.0),
+        (0.0, 0.0, 18.0),
+        (10.0, 8.0, 0.0),
+        (0.0, 11.0, 9.0),
+        (7.0, 0.0, 13.0),
+    )
+    pairs: list[MotionPair] = []
+    for session_index, (session_id, bias0) in enumerate(bias0_by_session.items()):
+        for pair_index, vector_deg in enumerate(vectors_deg):
+            r_b = so3_exp(np.deg2rad(np.asarray(vector_deg)))
+            r_a = r_true @ r_b @ r_true.T
+            index = session_index * 100 + pair_index
+            pairs.append(
+                MotionPair(
+                    session_id=session_id,
+                    i=index,
+                    j=index + 1,
+                    t_i_s=float(pair_index),
+                    t_j_s=float(pair_index + 1),
+                    R_A=r_a,
+                    R_B=r_b,
+                    t_A_m=np.zeros(3),
+                    t_B_m=np.zeros(3),
+                    metadata={
+                        "J_bg": np.eye(3).tolist(),
+                        "cov": (np.eye(3) * 1e-4).tolist(),
+                        "gyro_bias0_rad_s": bias0.tolist(),
+                    },
+                )
+            )
+
+    result = solve_joint_extrinsic(
+        pairs,
+        r_true,
+        force_rotation_only=True,
+        gyro_bias_rad_s_by_session=bias0_by_session,
+    )
+
+    assert result.phase_a_accepted
+    assert set(result.phase_a_comparison["variants"]) == {
+        "A0_fixed_bg_data_only",
+        "A1_session_bg_data_only",
+        "A2_session_bg_with_rotation_prior",
+    }
+    assert set(result.gyro_bias_rad_s_per_session) == {"s0", "s1"}
+    for session_id, bias0 in bias0_by_session.items():
+        np.testing.assert_allclose(
+            result.gyro_bias_rad_s_per_session[session_id], bias0, atol=1e-8
+        )
+    assert np.linalg.norm(so3_log(r_true.T @ result.T_IMU_lidar[:3, :3])) < 1e-8

tests/test_v1_pipeline.py

@@ -54,7 +54,7 @@ def _log(rotation: np.ndarray) -> np.ndarray:
     return so3_log(rotation)
 
 
-def test_synthetic_pipeline_rotation_and_time_offset(tmp_path: Path):
+def test_synthetic_pipeline_rejects_noisy_icp_but_keeps_time_audit(tmp_path: Path):
     meta = generate_synthetic_session(tmp_path, delta_t_s=0.17, yaw_extrinsic_deg=25.0)
     config = Path(__file__).resolve().parents[1] / "config" / "vehicle_installation.template.yaml"
     out = tmp_path / "out"
@@ -74,17 +74,28 @@ def test_synthetic_pipeline_rotation_and_time_offset(tmp_path: Path):
         min_pair_rotation_deg=2.0,
         min_pair_translation_m=0.05,
     )
-    result = run_calibration(request)
-    assert result.status.value == "rotation_only_accepted"
-    assert result.time_offset_s is not None
-    assert abs(result.time_offset_s - meta["delta_t_s"]) < 0.05
-    assert result.T_IMU_lidar is not None
-    # End-to-end uses approximate ICP; allow moderate absolute error but require consistency.
-    r_true = so3_exp(np.deg2rad(np.array([2.0, -1.5, meta["yaw_extrinsic_deg"]])))
-    err_deg = np.degrees(np.linalg.norm(_log(r_true.T @ result.T_IMU_lidar[:3, :3])))
-    assert err_deg < 15.0
+    progress_events: list[dict] = []
+    result = run_calibration(request, progress_callback=progress_events.append)
+    # The lightweight synthetic point cloud uses approximate ICP and has a
+    # roughly 3-degree P95 residual. The production gate must reject it rather
+    # than expose a plausible-looking extrinsic.
+    assert result.status.value == "blocked"
+    assert result.T_IMU_lidar is None
     session0 = result.details["sessions"][0]
-    assert session0["handeye"]["residual_rms_deg"] < 5.0
+    assert abs(session0["time_offset_s"] - meta["delta_t_s"]) < 0.05
+    assert result.details["joint_handeye"]["residual_p95_deg"] > 1.5
+    assert not result.details["joint_handeye"]["ok"]
+    assert progress_events[0]["event"] == "pipeline_start"
+    assert any(
+        event["stage"] == "motion_pairs" and event["event"] == "complete"
+        for event in progress_events
+    )
+    assert any(
+        event["stage"] == "joint_optimizer" and event["event"] == "phase_a_complete"
+        for event in progress_events
+    )
+    assert progress_events[-1]["stage"] == "finalize"
+    assert progress_events[-1]["event"] == "complete"
 
 
 def test_time_offset_on_synthetic(tmp_path: Path):
@@ -176,13 +187,17 @@ def test_synthetic_pipeline_full_se3_smoke(tmp_path: Path):
         "full_se3_accepted",
         "full_se3_rejected_due_to_observability",
         "rotation_only_accepted",
+        "blocked",
     }
-    assert result.T_IMU_lidar is not None
     session0 = result.details["sessions"][0]
     assert "delta_v" in session0.get("pair_notes", []) or session0.get("pair_count", 0) >= 0
-    # Phase-C fields appear when joint ran successfully on pairs.
-    if session0.get("ok"):
+    # Phase-C fields appear only when the strict rotation gate passed.
+    if result.status.value != "blocked":
+        assert result.T_IMU_lidar is not None
         assert "gyro_bias_rad_s" in session0["joint"]
+    else:
+        assert result.T_IMU_lidar is None
+        assert not result.details["joint_handeye"]["ok"]
 
 
 def test_signed_time_offset_refine_improves_or_keeps(tmp_path: Path):

tools/export_motion_pairs_for_viz.py

@@ -0,0 +1,111 @@
+#!/usr/bin/env python3
+"""Build motion_pairs.json next to an existing summary without re-solving extrinsic.
+
+Use this once for older calibration outputs that predate automatic pair caching.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import sys
+from pathlib import Path
+
+import numpy as np
+
+ROOT = Path(__file__).resolve().parents[1]
+if str(ROOT) not in sys.path:
+    sys.path.insert(0, str(ROOT))
+
+from imu_lidar.imu_audit import audit_imu
+from imu_lidar.imu_io import load_imu_samples
+from imu_lidar.keyframes import build_keyframes
+from imu_lidar.lidar_io import load_lidar_frames
+from imu_lidar.motion_pairs import build_motion_pairs
+from imu_lidar.motion_pairs_io import build_motion_pairs_payload, save_motion_pairs
+
+
+def _load_summary_meta(summary_path: Path) -> tuple[float, np.ndarray, str]:
+    summary = json.loads(summary_path.read_text(encoding="utf-8"))
+    delta_t = float(summary.get("time_offset_s") or 0.0)
+    session = (summary.get("details") or {}).get("sessions", [{}])[0]
+    session_id = str(session.get("session_id") or summary_path.parent.name)
+    bias = np.asarray(
+        (session.get("imu_audit") or {}).get("gyro_bias_rad_s")
+        or (session.get("joint") or {}).get("gyro_bias_rad_s")
+        or [0.0, 0.0, 0.0],
+        dtype=float,
+    ).reshape(3)
+    return delta_t, bias, session_id
+
+
+def export_one(
+    *,
+    lidar: Path,
+    imu: Path,
+    summary: Path,
+    output: Path | None,
+    min_rotation_deg: float,
+    min_translation_m: float,
+) -> Path:
+    delta_t, bias_from_summary, session_id = _load_summary_meta(summary)
+    imu_series = load_imu_samples(imu)
+    # Prefer freshly audited bias if summary bias is missing/zeros.
+    if float(np.linalg.norm(bias_from_summary)) < 1e-12:
+        bias = audit_imu(imu_series).gyro_bias_rad_s
+    else:
+        bias = bias_from_summary
+
+    frames = load_lidar_frames(lidar)
+    keyframes = build_keyframes(
+        frames,
+        min_translation_m=min_translation_m,
+        min_rotation_deg=min_rotation_deg,
+    )
+    pair_set = build_motion_pairs(
+        session_id=session_id,
+        keyframes=list(keyframes.frames),
+        keyframe_indices=keyframes.indices,
+        imu=imu_series,
+        delta_t_s=delta_t,
+        gyro_bias_rad_s=bias,
+        min_rotation_deg=min_rotation_deg,
+        min_translation_m=min_translation_m,
+    )
+    prepared = [
+        {
+            "session_id": session_id,
+            "time_offset_s": delta_t,
+            "gyro_bias_rad_s": np.asarray(bias, dtype=float).reshape(3),
+            "pairs": pair_set.pairs,
+        }
+    ]
+    payload = build_motion_pairs_payload(prepared_sessions=prepared)
+    out = output or (summary.parent / "motion_pairs.json")
+    save_motion_pairs(out, payload)
+    print(f"wrote {out} ({len(pair_set.pairs)} pairs, session={session_id}, dt={delta_t:.6f})")
+    return out
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--lidar", type=Path, required=True)
+    parser.add_argument("--imu", type=Path, required=True)
+    parser.add_argument("--summary", type=Path, required=True)
+    parser.add_argument("--output", type=Path, default=None, help="Default: <summary_dir>/motion_pairs.json")
+    parser.add_argument("--min-pair-rotation-deg", type=float, default=2.0)
+    parser.add_argument("--min-pair-translation-m", type=float, default=0.3)
+    args = parser.parse_args()
+    export_one(
+        lidar=args.lidar,
+        imu=args.imu,
+        summary=args.summary,
+        output=args.output,
+        min_rotation_deg=args.min_pair_rotation_deg,
+        min_translation_m=args.min_pair_translation_m,
+    )
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

tools/export_rscap_to_v1.py

@@ -1,13 +1,18 @@
 #!/usr/bin/env python3
-"""Export N300 IMU + H32 LiDAR captures to Lidar-IMU V1 intermediate format.
+"""Export IMU + H32 LiDAR captures to Lidar-IMU V1 intermediate format.
 
-Supported LiDAR sources (exactly one required):
+IMU sources:
+- ``--imu-kind hi13`` (HI13R4 / HI91) or ``n300`` or ``auto``
+- one or more ``--imu-rscap`` files (concatenated)
 
-- ``--lidar-dlog``: Medulla dlog from ``RSLidarH32_3D_DLogCaptureNet48``
-  (raw MSOP + DIFOP DObjects; preferred for new recordings)
-- ``--lidar-rscap``: legacy H32 MSOP V2 ``.rscap`` (MSOP-only defaults for angles)
+LiDAR sources (exactly one):
+- ``--lidar-dlog``: Medulla dlog dir **or recovered zip** (MSOP + DIFOP)
+- ``--lidar-rscap``: legacy H32 MSOP V2 ``.rscap``
 
-Output layout under --out:
+Optional host-time window (local wall clock, DateTime.Now.Ticks convention):
+- ``--host-start`` / ``--host-end`` e.g. ``2026-08-08T17:40:05``
+
+Output under ``--out``:
 
     imu.csv
     lidar/
@@ -15,8 +20,9 @@ Output layout under --out:
       frames/frame_XXXXX.npz
     export_summary.json
 
-Timestamps written into the intermediate format are **device times**
-(N300 device_timestamp_us, H32 MSOP device timestamp), not host receive time.
+Device times stay in ``t`` / ``t_start``/``t_end``. Host UTC receive times are
+also written so LiDAR–IMU alignment can bridge clocks without forcing first-frame
+device coincidence.
 """
 
 from __future__ import annotations
@@ -34,26 +40,48 @@ if str(ROOT) not in sys.path:
     sys.path.insert(0, str(ROOT))
 
 from tools.h32_dlog.load_session import load_h32_dlog_lidar
+from tools.h32_dlog.timeutil import (
+    local_wall_to_dotnet_ticks,
+    local_wall_to_utc_dotnet_ticks,
+    utc_dotnet_ticks_to_unix_s,
+)
 from tools.rscap_v2.capture_format_v2 import file_summary, read_capture
 from tools.rscap_v2.h32_msop import iter_h32_frames, iter_h32_frames_from_packets
-from tools.rscap_v2.n300_imu import iter_n300_imu_samples, samples_to_arrays
+from tools.rscap_v2.hi13_imu import iter_hi13_imu_samples
+from tools.rscap_v2.n300_imu import ImuSample, iter_n300_imu_samples, samples_to_arrays
 
 
-def write_imu_csv(path: Path, t: np.ndarray, gyro: np.ndarray, accel: np.ndarray) -> None:
+def write_imu_csv(path: Path, samples: list[ImuSample]) -> None:
     path.parent.mkdir(parents=True, exist_ok=True)
     with path.open("w", newline="", encoding="utf-8") as handle:
         writer = csv.writer(handle)
-        writer.writerow(["t", "gx", "gy", "gz", "ax", "ay", "az"])
-        for index in range(t.shape[0]):
+        writer.writerow(
+            [
+                "t",
+                "gx",
+                "gy",
+                "gz",
+                "ax",
+                "ay",
+                "az",
+                "t_host_utc_s",
+                "receive_utc_ticks",
+            ]
+        )
+        for sample in samples:
+            ticks = int(sample.host_receive_utc_ticks)
+            t_host = utc_dotnet_ticks_to_unix_s(ticks) if ticks > 0 else float("nan")
             writer.writerow(
                 [
-                    f"{t[index]:.9f}",
-                    f"{gyro[index, 0]:.12g}",
-                    f"{gyro[index, 1]:.12g}",
-                    f"{gyro[index, 2]:.12g}",
-                    f"{accel[index, 0]:.12g}",
-                    f"{accel[index, 1]:.12g}",
-                    f"{accel[index, 2]:.12g}",
+                    f"{sample.t_s:.9f}",
+                    f"{sample.gyro_rad_s[0]:.12g}",
+                    f"{sample.gyro_rad_s[1]:.12g}",
+                    f"{sample.gyro_rad_s[2]:.12g}",
+                    f"{sample.accel_m_s2[0]:.12g}",
+                    f"{sample.accel_m_s2[1]:.12g}",
+                    f"{sample.accel_m_s2[2]:.12g}",
+                    f"{t_host:.9f}" if ticks > 0 else "",
+                    ticks,
                 ]
             )
 
@@ -64,22 +92,45 @@ def write_lidar_session(root: Path, frames) -> dict:
     index_path = root / "frames_index.csv"
     with index_path.open("w", newline="", encoding="utf-8") as handle:
         writer = csv.writer(handle)
-        writer.writerow(["frame_id", "filename", "t_start", "t_end"])
+        writer.writerow(
+            [
+                "frame_id",
+                "filename",
+                "t_start",
+                "t_end",
+                "host_receive_utc_ticks",
+                "t_host_utc_s",
+                "host_receive_utc_end_ticks",
+                "t_host_utc_end_s",
+            ]
+        )
         point_counts = []
+        host_ok = 0
         for index, frame in enumerate(frames):
             rel = f"frames/frame_{index:05d}.npz"
             np.savez_compressed(root / rel, points=np.asarray(frame.points_xyz, dtype=np.float32))
+            h0 = int(getattr(frame, "host_receive_utc_ticks_start", 0) or 0)
+            h1 = int(getattr(frame, "host_receive_utc_ticks_end", 0) or 0)
+            t_host0 = utc_dotnet_ticks_to_unix_s(h0) if h0 > 0 else float("nan")
+            t_host1 = utc_dotnet_ticks_to_unix_s(h1) if h1 > 0 else float("nan")
+            if h0 > 0:
+                host_ok += 1
             writer.writerow(
                 [
                     index,
                     rel,
                     f"{frame.t_start_s:.9f}",
                     f"{frame.t_end_s:.9f}",
+                    h0,
+                    f"{t_host0:.9f}" if h0 > 0 else "",
+                    h1,
+                    f"{t_host1:.9f}" if h1 > 0 else "",
                 ]
             )
             point_counts.append(int(frame.points_xyz.shape[0]))
     return {
         "frames": len(frames),
+        "frames_with_host_utc": host_ok,
         "points_min": int(min(point_counts)) if point_counts else 0,
         "points_max": int(max(point_counts)) if point_counts else 0,
         "points_mean": float(np.mean(point_counts)) if point_counts else 0.0,
@@ -88,15 +139,63 @@ def write_lidar_session(root: Path, frames) -> dict:
     }
 
 
+def detect_imu_kind(paths: list[Path], explicit: str) -> str:
+    if explicit != "auto":
+        return explicit
+    joined = " ".join(path.name.lower() for path in paths)
+    if "hi13" in joined or "hipnuc" in joined:
+        return "hi13"
+    if "n300" in joined or "wheeltec" in joined:
+        return "n300"
+    return "hi13"
+
+
+def load_imu_samples(
+    paths: list[Path],
+    *,
+    kind: str,
+    host_ticks_min: int | None,
+    host_ticks_max: int | None,
+) -> tuple[list[ImuSample], list[dict], str]:
+    samples: list[ImuSample] = []
+    captures_meta: list[dict] = []
+    for path in paths:
+        capture = read_capture(path)
+        captures_meta.append(file_summary(capture))
+        if kind == "hi13":
+            part = iter_hi13_imu_samples(
+                capture,
+                host_utc_ticks_min=host_ticks_min,
+                host_utc_ticks_max=host_ticks_max,
+            )
+        elif kind == "n300":
+            part = iter_n300_imu_samples(capture)
+            if host_ticks_min is not None or host_ticks_max is not None:
+                part = [
+                    sample
+                    for sample in part
+                    if (host_ticks_min is None or sample.host_receive_utc_ticks >= host_ticks_min)
+                    and (host_ticks_max is None or sample.host_receive_utc_ticks <= host_ticks_max)
+                ]
+        else:
+            raise ValueError(f"unsupported imu kind: {kind}")
+        samples.extend(part)
+    samples.sort(key=lambda sample: (sample.t_s, sample.device_timestamp_us))
+    return samples, captures_meta, kind
+
+
 def export_session(
     *,
-    imu_rscap: Path,
+    imu_rscap: list[Path] | Path,
     out: Path,
     lidar_rscap: Path | None = None,
     lidar_dlog: Path | None = None,
+    imu_kind: str = "auto",
     msop_object: str = "frontlidar-msop-raw",
     difop_object: str = "frontlidar-difop-raw",
     require_difop: bool = False,
+    host_start: str | None = None,
+    host_end: str | None = None,
     frame_stride: int = 1,
     max_points_per_frame: int | None = 80000,
     min_range_m: float = 0.3,
@@ -106,24 +205,41 @@ def export_session(
     if (lidar_rscap is None) == (lidar_dlog is None):
         raise ValueError("provide exactly one of lidar_rscap or lidar_dlog")
 
+    imu_paths = [imu_rscap] if isinstance(imu_rscap, Path) else list(imu_rscap)
+    if not imu_paths:
+        raise ValueError("at least one --imu-rscap is required")
+
+    # LiDAR DObject tic uses DateTime.Now; IMU/MSOP host fields use UTC.
+    lidar_ticks_min = local_wall_to_dotnet_ticks(host_start) if host_start else None
+    lidar_ticks_max = local_wall_to_dotnet_ticks(host_end) if host_end else None
+    imu_ticks_min = local_wall_to_utc_dotnet_ticks(host_start) if host_start else None
+    imu_ticks_max = local_wall_to_utc_dotnet_ticks(host_end) if host_end else None
+    kind = detect_imu_kind(imu_paths, imu_kind)
+
     out.mkdir(parents=True, exist_ok=True)
-    imu_capture = read_capture(imu_rscap)
-
-    samples = iter_n300_imu_samples(imu_capture)
-    t, gyro, accel = samples_to_arrays(samples)
+    samples, imu_captures, kind = load_imu_samples(
+        imu_paths,
+        kind=kind,
+        host_ticks_min=imu_ticks_min,
+        host_ticks_max=imu_ticks_max,
+    )
+    t, _gyro, _accel = samples_to_arrays(samples)
     imu_csv = out / "imu.csv"
-    write_imu_csv(imu_csv, t, gyro, accel)
+    write_imu_csv(imu_csv, samples)
+    imu_host_ok = sum(1 for sample in samples if sample.host_receive_utc_ticks > 0)
 
-    lidar_meta: dict
     if lidar_dlog is not None:
         session = load_h32_dlog_lidar(
             lidar_dlog,
             msop_object=msop_object,
             difop_object=difop_object,
             require_difop=require_difop,
+            host_ticks_min=lidar_ticks_min,
+            host_ticks_max=lidar_ticks_max,
         )
         frames = iter_h32_frames_from_packets(
             session.msop_packets,
+            host_utc_ticks=session.msop_host_utc_ticks,
             min_frame_points=min_frame_points,
             frame_stride=frame_stride,
             min_range_m=min_range_m,
@@ -134,16 +250,20 @@ def export_session(
         )
         lidar_meta = {
             "source": "dlog",
-            "lidar_dlog": str(session.dlog_root),
+            "lidar_dlog": session.dlog_root,
             "msop_object": session.msop_object,
             "difop_object": session.difop_object,
             "msop_packets": len(session.msop_packets),
+            "msop_packets_with_host_utc": sum(1 for ticks in session.msop_host_utc_ticks if ticks > 0),
             "msop_batches": session.msop_batch_count,
             "difop_records": session.difop_record_count,
             "session_id": session.session_id,
             "lidar_ip": session.lidar_ip,
             "angle_source": session.angle_source,
-            "timestamp_note": "h32_msop_device_timestamp -> seconds (from MSOP bytes)",
+            "timestamp_note": (
+                "device: h32_msop_device_timestamp -> seconds; "
+                "host: MSOP HostReceiveUtcTicks -> unix seconds"
+            ),
         }
     else:
         assert lidar_rscap is not None
@@ -161,25 +281,46 @@ def export_session(
             "lidar_rscap": str(lidar_rscap),
             "capture": file_summary(lidar_capture),
             "angle_source": "default_msop_only_vertical_-16_to_16_deg",
-            "timestamp_note": "h32_msop_device_timestamp_ms -> seconds",
+            "timestamp_note": (
+                "device: h32_msop_device_timestamp_ms -> seconds; "
+                "host: rscap receive_utc_ticks -> unix seconds"
+            ),
         }
 
     lidar_dir = out / "lidar"
     lidar_stats = write_lidar_session(lidar_dir, frames)
+    imu_time_note = (
+        "hi13_device_timestamp_ms -> seconds"
+        if kind == "hi13"
+        else "n300_device_timestamp_us -> seconds"
+    )
 
     summary = {
-        "imu_rscap": str(imu_rscap),
+        "imu_rscap": [str(path) for path in imu_paths],
+        "imu_kind": kind,
         "out": str(out),
+        "host_window": {
+            "host_start": host_start,
+            "host_end": host_end,
+            "lidar_ticks_min": lidar_ticks_min,
+            "lidar_ticks_max": lidar_ticks_max,
+            "imu_ticks_min": imu_ticks_min,
+            "imu_ticks_max": imu_ticks_max,
+            "note": "local wall cut; lidar DObject tic=DateTime.Now, IMU/MSOP host=UTC",
+        },
         "timestamp_policy": {
-            "imu": "n300_device_timestamp_us -> seconds",
-            "lidar": lidar_meta["timestamp_note"],
-            "host_utc": "not used as calibration timeline",
+            "imu_device": imu_time_note,
+            "imu_host": "rscap receive_utc_ticks -> t_host_utc_s",
+            "lidar_device": "MSOP device timestamp -> t_start/t_end",
+            "lidar_host": "MSOP HostReceiveUtcTicks -> t_host_utc_s",
+            "calibration_align": "bridge via host UTC; do not force first device samples to coincide",
         },
         "imu": {
             "samples": int(t.shape[0]),
+            "samples_with_host_utc": imu_host_ok,
             "t_start": float(t[0]) if t.size else None,
             "t_end": float(t[-1]) if t.size else None,
-            "capture": file_summary(imu_capture),
+            "captures": imu_captures,
         },
         "lidar": {
             **lidar_stats,
@@ -201,41 +342,38 @@ def export_session(
 
 def main() -> int:
     parser = argparse.ArgumentParser(description=__doc__)
-    parser.add_argument("--imu-rscap", type=Path, required=True, help="N300 V2 .rscap")
+    parser.add_argument(
+        "--imu-rscap",
+        type=Path,
+        action="append",
+        required=True,
+        help="IMU V2 .rscap (repeatable)",
+    )
+    parser.add_argument(
+        "--imu-kind",
+        choices=("auto", "hi13", "n300"),
+        default="auto",
+        help="IMU decoder (default: auto from filename)",
+    )
     lidar = parser.add_mutually_exclusive_group(required=True)
     lidar.add_argument(
         "--lidar-dlog",
         type=Path,
-        help="H32 Medulla dlog root (dobject/ + dobject_recording/), preferred",
+        help="H32 dlog directory or recovered zip (indices.log + data.bin)",
     )
     lidar.add_argument(
         "--lidar-rscap",
         type=Path,
-        help="Legacy H32 MSOP V2 .rscap (no DIFOP; default vertical angles)",
-    )
-    parser.add_argument(
-        "--msop-object",
-        default="frontlidar-msop-raw",
-        help="DObject name for raw MSOP batches (dlog path)",
-    )
-    parser.add_argument(
-        "--difop-object",
-        default="frontlidar-difop-raw",
-        help="DObject name for raw DIFOP packets (dlog path)",
-    )
-    parser.add_argument(
-        "--require-difop",
-        action="store_true",
-        help="Fail if dlog has no valid DIFOP channel angles",
-    )
-    parser.add_argument("--out", type=Path, required=True, help="Output session directory")
-    parser.add_argument("--frame-stride", type=int, default=1, help="Keep every N-th LiDAR frame")
-    parser.add_argument(
-        "--max-points-per-frame",
-        type=int,
-        default=80000,
-        help="Uniform downsample cap per frame; 0 disables",
+        help="Legacy H32 MSOP V2 .rscap",
     )
+    parser.add_argument("--msop-object", default="frontlidar-msop-raw")
+    parser.add_argument("--difop-object", default="frontlidar-difop-raw")
+    parser.add_argument("--require-difop", action="store_true")
+    parser.add_argument("--host-start", type=str, default=None, help="Local wall start, e.g. 2026-08-08T17:40:05")
+    parser.add_argument("--host-end", type=str, default=None, help="Local wall end, e.g. 2026-08-08T17:45:15")
+    parser.add_argument("--out", type=Path, required=True)
+    parser.add_argument("--frame-stride", type=int, default=1)
+    parser.add_argument("--max-points-per-frame", type=int, default=80000)
     parser.add_argument("--min-range-m", type=float, default=0.3)
     parser.add_argument("--max-range-m", type=float, default=120.0)
     parser.add_argument("--min-frame-points", type=int, default=100)
@@ -245,9 +383,12 @@ def main() -> int:
         imu_rscap=args.imu_rscap,
         lidar_rscap=args.lidar_rscap,
         lidar_dlog=args.lidar_dlog,
+        imu_kind=args.imu_kind,
         msop_object=args.msop_object,
         difop_object=args.difop_object,
         require_difop=args.require_difop,
+        host_start=args.host_start,
+        host_end=args.host_end,
         out=args.out,
         frame_stride=args.frame_stride,
         max_points_per_frame=max_points,
@@ -258,10 +399,14 @@ def main() -> int:
     print(
         json.dumps(
             {
+                "imu_kind": summary["imu_kind"],
                 "imu_samples": summary["imu"]["samples"],
+                "imu_host_utc": summary["imu"]["samples_with_host_utc"],
                 "lidar_frames": summary["lidar"]["frames"],
+                "lidar_host_utc": summary["lidar"]["frames_with_host_utc"],
                 "lidar_source": summary["lidar"]["source"],
                 "angle_source": summary["lidar"]["angle_source"],
+                "host_window": summary["host_window"],
                 "imu_csv": summary["outputs"]["imu_csv"],
                 "lidar_session": summary["outputs"]["lidar_session"],
                 "export_summary": str(Path(args.out) / "export_summary.json"),
@@ -271,9 +416,13 @@ def main() -> int:
         )
     )
     if summary["imu"]["samples"] == 0:
-        raise SystemExit("no valid N300 IMU samples decoded")
+        raise SystemExit("no valid IMU samples decoded in window")
     if summary["lidar"]["frames"] == 0:
-        raise SystemExit("no valid H32 frames decoded")
+        raise SystemExit("no valid H32 frames decoded in window")
+    if summary["lidar"]["frames_with_host_utc"] == 0:
+        raise SystemExit("no LiDAR frames with MSOP HostReceiveUtcTicks; cannot host-bridge align")
+    if summary["imu"]["samples_with_host_utc"] == 0:
+        raise SystemExit("no IMU samples with host receive UTC; cannot host-bridge align")
     return 0
 
 

tools/export_usable_20260808_windows.py

@@ -0,0 +1,125 @@
+#!/usr/bin/env python3
+"""Export priority LiDAR–IMU windows from calibration_usable_20260808.
+
+Does not push anything; writes local V1 sessions under --out-root.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import sys
+from pathlib import Path
+
+ROOT = Path(__file__).resolve().parents[1]
+if str(ROOT) not in sys.path:
+    sys.path.insert(0, str(ROOT))
+
+from tools.export_rscap_to_v1 import export_session
+
+DEFAULT_DATA = Path(r"D:\data\calibration_usable_20260808")
+LIDAR_ZIP = "lidar_dlog/dorec_recovered_20260808_171438_181422.zip"
+IMU_MAIN = "imu_rscap/hi13r4-imu_20260808-092827.638_39783edb-e46e-4b28-a5e8-427b981c2fce.rscap"
+IMU_TAIL = "imu_rscap/hi13r4-imu_20260808-101022.036_87ea5edc-cd3d-4192-809a-469fbc8cac01.rscap"
+
+# From usable-segment chart (local wall clock).
+WINDOWS = [
+    {
+        "name": "priority_174005_174515",
+        "host_start": "2026-08-08T17:40:05",
+        "host_end": "2026-08-08T17:45:15",
+        "imu": [IMU_MAIN],
+        "priority": True,
+    },
+    {
+        "name": "priority_174905_175450",
+        "host_start": "2026-08-08T17:49:05",
+        "host_end": "2026-08-08T17:54:50",
+        "imu": [IMU_MAIN],
+        "priority": True,
+    },
+    {
+        "name": "priority_175910_180530",
+        "host_start": "2026-08-08T17:59:10",
+        "host_end": "2026-08-08T18:05:30",
+        "imu": [IMU_MAIN],
+        "priority": True,
+    },
+    {
+        "name": "usable_181035_181050",
+        "host_start": "2026-08-08T18:10:35",
+        "host_end": "2026-08-08T18:10:50",
+        "imu": [IMU_TAIL],
+        "priority": False,
+    },
+    {
+        "name": "usable_181225_181300",
+        "host_start": "2026-08-08T18:12:25",
+        "host_end": "2026-08-08T18:13:00",
+        "imu": [IMU_TAIL],
+        "priority": False,
+    },
+    {
+        "name": "usable_181350_181410",
+        "host_start": "2026-08-08T18:13:50",
+        "host_end": "2026-08-08T18:14:10",
+        "imu": [IMU_TAIL],
+        "priority": False,
+    },
+]
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--data-root", type=Path, default=DEFAULT_DATA)
+    parser.add_argument(
+        "--out-root",
+        type=Path,
+        default=DEFAULT_DATA / "sessions_v1",
+    )
+    parser.add_argument("--priority-only", action="store_true", default=True)
+    parser.add_argument("--all-windows", action="store_true")
+    parser.add_argument("--frame-stride", type=int, default=5)
+    parser.add_argument("--max-points-per-frame", type=int, default=40000)
+    args = parser.parse_args()
+    priority_only = not args.all_windows
+    lidar = args.data_root / LIDAR_ZIP
+    if not lidar.is_file():
+        raise SystemExit(f"missing lidar zip: {lidar}")
+
+    selected = [w for w in WINDOWS if (not priority_only) or w["priority"]]
+    results = []
+    for window in selected:
+        out = args.out_root / window["name"]
+        imu_paths = [args.data_root / rel for rel in window["imu"]]
+        print(f"=== exporting {window['name']} ===", flush=True)
+        summary = export_session(
+            imu_rscap=imu_paths,
+            lidar_dlog=lidar,
+            imu_kind="hi13",
+            require_difop=True,
+            host_start=window["host_start"],
+            host_end=window["host_end"],
+            out=out,
+            frame_stride=args.frame_stride,
+            max_points_per_frame=args.max_points_per_frame,
+        )
+        brief = {
+            "name": window["name"],
+            "imu_samples": summary["imu"]["samples"],
+            "lidar_frames": summary["lidar"]["frames"],
+            "angle_source": summary["lidar"]["angle_source"],
+            "out": str(out),
+        }
+        results.append(brief)
+        print(json.dumps(brief, ensure_ascii=False, indent=2), flush=True)
+
+    manifest = args.out_root / "export_windows_manifest.json"
+    args.out_root.mkdir(parents=True, exist_ok=True)
+    manifest.write_text(json.dumps(results, ensure_ascii=False, indent=2) + "\n", encoding="utf-8")
+    print(f"manifest: {manifest}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

tools/h32_dlog/__init__.py

@@ -1,12 +1,18 @@
 """Medulla dlog readers for RSLidarH32_3D_DLogCaptureNet48 raw MSOP/DIFOP."""
 
 from .difop import parse_difop_angles
-from .dobject import discover_records, iter_payloads, resolve_dlog_root
+from .dobject import discover_records, iter_payloads, open_dlog_source, resolve_dlog_root
+from .load_session import H32DlogLidarSession, load_h32_dlog_lidar
 from .payload_v1 import parse_difop_payload, parse_msop_batch_payload
+from .timeutil import local_wall_to_dotnet_ticks
 
 __all__ = [
+    "H32DlogLidarSession",
     "discover_records",
     "iter_payloads",
+    "load_h32_dlog_lidar",
+    "local_wall_to_dotnet_ticks",
+    "open_dlog_source",
     "parse_difop_angles",
     "parse_difop_payload",
     "parse_msop_batch_payload",

tools/h32_dlog/dobject.py

@@ -1,16 +1,24 @@
-"""Index and read Medulla DObject recordings (dobject/ + dobject_recording/)."""
+"""Index and read Medulla DObject recordings.
+
+Supports:
+
+- standard layout: ``dobject/**/*.log`` + ``dobject_recording/**/*.dorec``
+- recovered layout: ``dobject/all/indices.log`` + ``dobject_recording/data.bin``
+- either as an extracted directory or a zip containing those paths
+"""
 
 from __future__ import annotations
 
 import re
 import struct
+import zipfile
 from dataclasses import dataclass
 from pathlib import Path
 from typing import BinaryIO, Iterator
 
 
 RECORD_RE = re.compile(
-    r"^\[(?P<log_time>[^]]+)\].*?DObject `(?P<name>[^`]+)` post "
+    r"^(?:\[(?P<log_time>[^]]+)\])?>?\s*DObject `(?P<name>[^`]+)` post "
     r"len=(?P<len>\d+)B, id:(?P<id>[0-9A-Fa-f]+), tic:(?P<tic>\d+), "
     r"@(?P<file>[^:]+):(?P<offset>\d+)"
 )
@@ -29,37 +37,225 @@ class RecordRef:
     dotnet_ticks: int
 
 
-def resolve_dlog_root(value: Path | str) -> Path:
-    root = Path(value).expanduser().resolve()
-    if (root / "dobject").is_dir() and (root / "dobject_recording").is_dir():
-        return root
-    child = root / "dlog"
-    if (child / "dobject").is_dir() and (child / "dobject_recording").is_dir():
-        return child
-    raise FileNotFoundError(f"{root} does not contain dobject and dobject_recording")
+class _ZipStoredMemberIO:
+    """Random-access reader for a ZIP_STORED member via the underlying zip file.
+
+    ``ZipExtFile.seek`` on multi-GB members is far too slow for per-record reads.
+    """
+
+    def __init__(self, zip_path: Path, member_name: str, data_offset: int, data_size: int):
+        self._path = zip_path
+        self._member_name = member_name
+        self._data_offset = data_offset
+        self._data_size = data_size
+        self._fh = zip_path.open("rb")
+        self._pos = 0
+
+    def seek(self, offset: int, whence: int = 0) -> int:
+        if whence == 0:
+            self._pos = offset
+        elif whence == 1:
+            self._pos += offset
+        elif whence == 2:
+            self._pos = self._data_size + offset
+        else:
+            raise ValueError(f"invalid whence: {whence}")
+        if self._pos < 0:
+            raise ValueError("negative seek")
+        return self._pos
+
+    def read(self, size: int = -1) -> bytes:
+        if size is None or size < 0:
+            size = self._data_size - self._pos
+        if size <= 0 or self._pos >= self._data_size:
+            return b""
+        size = min(size, self._data_size - self._pos)
+        self._fh.seek(self._data_offset + self._pos)
+        data = self._fh.read(size)
+        self._pos += len(data)
+        return data
+
+    def close(self) -> None:
+        self._fh.close()
 
 
-def discover_records(dlog_root: Path, object_name: str) -> list[RecordRef]:
-    pending: list[tuple[str, str, str, int, int, str, int, str]] = []
-    for log_path in sorted((dlog_root / "dobject").rglob("*.log")):
-        relative_log = log_path.relative_to(dlog_root).as_posix()
-        with log_path.open("r", encoding="utf-8", errors="replace") as stream:
-            for line in stream:
-                match = RECORD_RE.search(line)
-                if not match or match.group("name").casefold() != object_name.casefold():
-                    continue
-                pending.append(
-                    (
-                        match.group("name"),
-                        match.group("log_time"),
-                        relative_log,
-                        int(match.group("offset")),
-                        int(match.group("len")),
-                        match.group("id").upper(),
-                        int(match.group("tic")),
-                        match.group("file"),
-                    )
+def _zip_stored_member_offset(zip_path: Path, info: zipfile.ZipInfo) -> int:
+    if info.compress_type != zipfile.ZIP_STORED:
+        raise RuntimeError(
+            f"member {info.filename!r} is compressed (type={info.compress_type}); "
+            "extract it first or store uncompressed"
+        )
+    with zip_path.open("rb") as handle:
+        handle.seek(info.header_offset)
+        header = handle.read(30)
+        if len(header) != 30 or header[:4] != b"PK\x03\x04":
+            raise RuntimeError(f"bad local zip header for {info.filename!r}")
+        name_len, extra_len = struct.unpack("<HH", header[26:30])
+        return info.header_offset + 30 + name_len + extra_len
+
+
+@dataclass
+class DlogSource:
+    """Opened dlog directory or recovered zip."""
+
+    label: str
+    directory: Path | None = None
+    zip_path: Path | None = None
+    _zip: zipfile.ZipFile | None = None
+    _log_cache: dict[str, str] | None = None
+    _member_offsets: dict[str, tuple[int, int]] | None = None
+
+    def close(self) -> None:
+        if self._zip is not None:
+            self._zip.close()
+            self._zip = None
+
+    def __enter__(self) -> "DlogSource":
+        return self
+
+    def __exit__(self, exc_type, exc, tb) -> None:
+        self.close()
+
+    def iter_log_texts(self) -> Iterator[tuple[str, str]]:
+        if self.zip_path is not None:
+            assert self._zip is not None
+            if self._log_cache is None:
+                self._log_cache = {}
+                names = sorted(
+                    name
+                    for name in self._zip.namelist()
+                    if name.replace("\\", "/").startswith("dobject/")
+                    and name.replace("\\", "/").endswith(".log")
                 )
+                for name in names:
+                    key = name.replace("\\", "/")
+                    self._log_cache[key] = self._zip.read(name).decode("utf-8", errors="replace")
+            for name, text in self._log_cache.items():
+                yield name, text
+            return
+        assert self.directory is not None
+        for log_path in sorted((self.directory / "dobject").rglob("*.log")):
+            relative = log_path.relative_to(self.directory).as_posix()
+            yield relative, log_path.read_text(encoding="utf-8", errors="replace")
+
+    def open_recording(self, name: str) -> tuple[object, BinaryIO]:
+        """Return (owner, binary stream) supporting seek/read of one recording member."""
+
+        base = Path(name).name
+        if self.zip_path is not None:
+            assert self._zip is not None
+            candidates = [
+                n
+                for n in self._zip.namelist()
+                if Path(n.replace("\\", "/")).name.casefold() == base.casefold()
+                and "dobject_recording/" in n.replace("\\", "/")
+            ]
+            if not candidates:
+                alt = name.replace("\\", "/")
+                if alt in self._zip.namelist():
+                    candidates = [alt]
+                elif f"dobject_recording/{base}" in self._zip.namelist():
+                    candidates = [f"dobject_recording/{base}"]
+            if not candidates:
+                raise FileNotFoundError(f"missing recording in zip: {name}")
+            if len(candidates) > 1:
+                raise RuntimeError(f"ambiguous recording in zip {name}: {candidates}")
+            member = candidates[0].replace("\\", "/")
+            if self._member_offsets is None:
+                self._member_offsets = {}
+            if member not in self._member_offsets:
+                info = self._zip.getinfo(member)
+                self._member_offsets[member] = (
+                    _zip_stored_member_offset(self.zip_path, info),
+                    info.file_size,
+                )
+            data_offset, data_size = self._member_offsets[member]
+            stream = _ZipStoredMemberIO(self.zip_path, member, data_offset, data_size)
+            return stream, stream
+
+        assert self.directory is not None
+        index = index_dorec_files(self.directory)
+        if base.casefold() == "data.bin":
+            path = self.directory / "dobject_recording" / "data.bin"
+            if not path.is_file():
+                matches = list((self.directory / "dobject_recording").rglob("data.bin"))
+                if not matches:
+                    raise FileNotFoundError(f"missing recording file: {name}")
+                path = matches[0]
+            stream = path.open("rb")
+            return stream, stream
+        path = choose_dorec(index, name)
+        stream = path.open("rb")
+        return stream, stream
+
+
+def open_dlog_source(value: Path | str) -> DlogSource:
+    path = Path(value).expanduser().resolve()
+    if path.is_file() and path.suffix.lower() == ".zip":
+        zf = zipfile.ZipFile(path, "r")
+        names = {n.replace("\\", "/") for n in zf.namelist()}
+        has_log = any(n.startswith("dobject/") and n.endswith(".log") for n in names)
+        has_rec = any(n.startswith("dobject_recording/") for n in names)
+        if not (has_log and has_rec):
+            zf.close()
+            raise FileNotFoundError(f"{path} is not a recovered/standard dlog zip")
+        return DlogSource(label=str(path), zip_path=path, _zip=zf)
+
+    root = path
+    if not ((root / "dobject").is_dir() and (root / "dobject_recording").is_dir()):
+        child = root / "dlog"
+        if (child / "dobject").is_dir() and (child / "dobject_recording").is_dir():
+            root = child
+        else:
+            raise FileNotFoundError(f"{path} does not contain dobject and dobject_recording")
+    return DlogSource(label=str(root), directory=root)
+
+
+def resolve_dlog_root(value: Path | str) -> Path:
+    """Backward-compatible helper: directory roots only (not zip)."""
+
+    source = open_dlog_source(value)
+    try:
+        if source.directory is None:
+            raise FileNotFoundError(
+                f"{value} is a zip; use open_dlog_source()/iter_payloads_from_source()"
+            )
+        return source.directory
+    finally:
+        source.close()
+
+
+def discover_records_from_source(
+    source: DlogSource,
+    object_name: str,
+    *,
+    host_ticks_min: int | None = None,
+    host_ticks_max: int | None = None,
+) -> list[RecordRef]:
+    pending: list[tuple[str, str, str, int, int, str, int, str]] = []
+    name_key = object_name.casefold()
+    for relative_log, text in source.iter_log_texts():
+        for line in text.splitlines():
+            match = RECORD_RE.search(line.strip())
+            if not match or match.group("name").casefold() != name_key:
+                continue
+            ticks = int(match.group("tic"))
+            if host_ticks_min is not None and ticks < host_ticks_min:
+                continue
+            if host_ticks_max is not None and ticks > host_ticks_max:
+                continue
+            pending.append(
+                (
+                    match.group("name"),
+                    match.group("log_time") or "",
+                    relative_log,
+                    int(match.group("offset")),
+                    int(match.group("len")),
+                    match.group("id").upper(),
+                    ticks,
+                    match.group("file"),
+                )
+            )
     pending.sort(key=lambda item: (item[6], item[7].casefold(), item[3]))
     seen: set[tuple[str, int, int]] = set()
     records: list[RecordRef] = []
@@ -84,10 +280,19 @@ def discover_records(dlog_root: Path, object_name: str) -> list[RecordRef]:
     return records
 
 
+def discover_records(dlog_root: Path, object_name: str) -> list[RecordRef]:
+    with open_dlog_source(dlog_root) as source:
+        return discover_records_from_source(source, object_name)
+
+
 def index_dorec_files(dlog_root: Path) -> dict[str, list[Path]]:
     result: dict[str, list[Path]] = {}
-    for path in (dlog_root / "dobject_recording").rglob("*.dorec"):
-        result.setdefault(path.name.casefold(), []).append(path)
+    recording = dlog_root / "dobject_recording"
+    if not recording.is_dir():
+        return result
+    for path in recording.rglob("*"):
+        if path.is_file() and path.suffix.lower() in {".dorec", ".bin"}:
+            result.setdefault(path.name.casefold(), []).append(path)
     return result
 
 
@@ -107,17 +312,15 @@ def read_exact(stream: BinaryIO, size: int) -> bytes:
     return data
 
 
-def read_record_payload(path: Path, record: RecordRef) -> bytes:
-    with path.open("rb") as stream:
-        stream.seek(record.source_offset)
-        name_length = read_exact(stream, 1)[0]
-        name = read_exact(stream, name_length).decode("ascii")
-        ticks = struct.unpack("<q", read_exact(stream, 8))[0]
-        id_length = read_exact(stream, 1)[0]
-        id_bytes = read_exact(stream, id_length)
-        payload_length = struct.unpack("<i", read_exact(stream, 4))[0]
-        payload = read_exact(stream, payload_length)
-
+def _read_payload_at(stream: BinaryIO, record: RecordRef) -> bytes:
+    stream.seek(record.source_offset)
+    name_length = read_exact(stream, 1)[0]
+    name = read_exact(stream, name_length).decode("ascii")
+    ticks = struct.unpack("<q", read_exact(stream, 8))[0]
+    id_length = read_exact(stream, 1)[0]
+    id_bytes = read_exact(stream, id_length)
+    payload_length = struct.unpack("<i", read_exact(stream, 4))[0]
+    payload = read_exact(stream, payload_length)
     try:
         record_id = id_bytes.decode("ascii")
     except UnicodeDecodeError:
@@ -133,47 +336,46 @@ def read_record_payload(path: Path, record: RecordRef) -> bytes:
     return payload
 
 
-def iter_payloads(dlog_root: Path, object_name: str) -> Iterator[tuple[RecordRef, bytes]]:
-    root = resolve_dlog_root(dlog_root)
-    records = discover_records(root, object_name)
+def iter_payloads_from_source(
+    source: DlogSource,
+    object_name: str,
+    *,
+    host_ticks_min: int | None = None,
+    host_ticks_max: int | None = None,
+) -> Iterator[tuple[RecordRef, bytes]]:
+    records = discover_records_from_source(
+        source,
+        object_name,
+        host_ticks_min=host_ticks_min,
+        host_ticks_max=host_ticks_max,
+    )
     if not records:
         return
-    dorec_index = index_dorec_files(root)
-    open_files: dict[str, tuple[Path, BinaryIO]] = {}
+    open_files: dict[str, BinaryIO] = {}
     try:
         for record in records:
-            key = record.source_dorec.casefold()
-            handle = open_files.get(key)
-            if handle is None:
-                path = choose_dorec(dorec_index, record.source_dorec)
-                handle = (path, path.open("rb"))
-                open_files[key] = handle
-            path, stream = handle
-            stream.seek(record.source_offset)
-            name_length = read_exact(stream, 1)[0]
-            name = read_exact(stream, name_length).decode("ascii")
-            ticks = struct.unpack("<q", read_exact(stream, 8))[0]
-            id_length = read_exact(stream, 1)[0]
-            id_bytes = read_exact(stream, id_length)
-            payload_length = struct.unpack("<i", read_exact(stream, 4))[0]
-            payload = read_exact(stream, payload_length)
-            try:
-                record_id = id_bytes.decode("ascii")
-            except UnicodeDecodeError:
-                record_id = id_bytes.hex().upper()
-            if name != record.object_name:
-                raise ValueError(f"name mismatch: log={record.object_name}, dorec={name}")
-            if ticks != record.dotnet_ticks:
-                raise ValueError(f"tick mismatch: log={record.dotnet_ticks}, dorec={ticks}")
-            if payload_length != record.payload_length:
-                raise ValueError(
-                    f"payload mismatch: log={record.payload_length}, dorec={payload_length}"
-                )
-            if record_id.upper() != record.log_record_id.upper():
-                raise ValueError(
-                    f"record id mismatch: log={record.log_record_id}, dorec={record_id}"
-                )
-            yield record, payload
+            key = Path(record.source_dorec).name.casefold()
+            stream = open_files.get(key)
+            if stream is None:
+                _owner, stream = source.open_recording(record.source_dorec)
+                open_files[key] = stream
+            yield record, _read_payload_at(stream, record)
     finally:
-        for _path, stream in open_files.values():
+        for stream in open_files.values():
             stream.close()
+
+
+def iter_payloads(
+    dlog_root: Path | str,
+    object_name: str,
+    *,
+    host_ticks_min: int | None = None,
+    host_ticks_max: int | None = None,
+) -> Iterator[tuple[RecordRef, bytes]]:
+    with open_dlog_source(dlog_root) as source:
+        yield from iter_payloads_from_source(
+            source,
+            object_name,
+            host_ticks_min=host_ticks_min,
+            host_ticks_max=host_ticks_max,
+        )

tools/h32_dlog/load_session.py

@@ -10,16 +10,21 @@ import numpy as np
 from tools.rscap_v2.h32_msop import default_horizontal_deg, default_vertical_deg
 
 from .difop import DifopAngles, parse_difop_angles
-from .dobject import discover_records, iter_payloads, resolve_dlog_root
+from .dobject import (
+    discover_records_from_source,
+    iter_payloads_from_source,
+    open_dlog_source,
+)
 from .payload_v1 import parse_difop_payload, parse_msop_batch_payload
 
 
 @dataclass
 class H32DlogLidarSession:
-    dlog_root: Path
+    dlog_root: str
     msop_object: str
     difop_object: str
     msop_packets: list[bytes]
+    msop_host_utc_ticks: list[int]
     msop_batch_count: int
     difop_record_count: int
     angle_source: str
@@ -27,6 +32,8 @@ class H32DlogLidarSession:
     horizontal_deg: np.ndarray
     session_id: str | None = None
     lidar_ip: str | None = None
+    host_ticks_min: int | None = None
+    host_ticks_max: int | None = None
 
 
 def load_h32_dlog_lidar(
@@ -35,65 +42,99 @@ def load_h32_dlog_lidar(
     msop_object: str = "frontlidar-msop-raw",
     difop_object: str = "frontlidar-difop-raw",
     require_difop: bool = False,
+    host_ticks_min: int | None = None,
+    host_ticks_max: int | None = None,
 ) -> H32DlogLidarSession:
-    root = resolve_dlog_root(dlog_root)
-    msop_packets: list[bytes] = []
-    batch_count = 0
-    session_id: str | None = None
-    lidar_ip: str | None = None
+    with open_dlog_source(dlog_root) as source:
+        # DIFOP angles: prefer packets inside the window, else any in the capture.
+        angles: DifopAngles | None = None
+        difop_count = 0
+        session_id: str | None = None
+        lidar_ip: str | None = None
+        for _record, payload in iter_payloads_from_source(
+            source,
+            difop_object,
+            host_ticks_min=host_ticks_min,
+            host_ticks_max=host_ticks_max,
+        ):
+            difop = parse_difop_payload(payload)
+            difop_count += 1
+            try:
+                angles = parse_difop_angles(difop.raw)
+            except ValueError:
+                continue
+            if session_id is None:
+                session_id = difop.session_id
+                lidar_ip = difop.lidar_ip
 
-    for _record, payload in iter_payloads(root, msop_object):
-        batch = parse_msop_batch_payload(payload)
-        batch_count += 1
-        if session_id is None:
-            session_id = batch.session_id
-            lidar_ip = batch.lidar_ip
-        for item in batch.packets:
-            msop_packets.append(item.raw)
+        if angles is None:
+            for _record, payload in iter_payloads_from_source(source, difop_object):
+                difop = parse_difop_payload(payload)
+                difop_count += 1
+                try:
+                    angles = parse_difop_angles(difop.raw)
+                except ValueError:
+                    continue
+                if session_id is None:
+                    session_id = difop.session_id
+                    lidar_ip = difop.lidar_ip
+                if angles is not None:
+                    break
 
-    angles: DifopAngles | None = None
-    difop_count = 0
-    for _record, payload in iter_payloads(root, difop_object):
-        difop = parse_difop_payload(payload)
-        difop_count += 1
-        try:
-            angles = parse_difop_angles(difop.raw)
-        except ValueError:
-            continue
-        if session_id is None:
-            session_id = difop.session_id
-            lidar_ip = difop.lidar_ip
+        msop_packets: list[bytes] = []
+        msop_host_utc_ticks: list[int] = []
+        batch_count = 0
+        for record, payload in iter_payloads_from_source(
+            source,
+            msop_object,
+            host_ticks_min=host_ticks_min,
+            host_ticks_max=host_ticks_max,
+        ):
+            batch = parse_msop_batch_payload(payload)
+            batch_count += 1
+            if session_id is None:
+                session_id = batch.session_id
+                lidar_ip = batch.lidar_ip
+            for item in batch.packets:
+                msop_packets.append(item.raw)
+                # Per-packet UTC host receive from MSOP DLog payload only.
+                # Do NOT fall back to DObject tic (DateTime.Now / local).
+                msop_host_utc_ticks.append(int(item.host_receive_utc_ticks))
 
-    if not msop_packets:
-        msop_records = discover_records(root, msop_object)
-        raise RuntimeError(
-            f"no MSOP packets from DObject {msop_object!r} under {root} "
-            f"(log records={len(msop_records)})"
-        )
-
-    if angles is None:
-        if require_difop:
+        if not msop_packets:
+            msop_records = discover_records_from_source(source, msop_object)
             raise RuntimeError(
-                f"no valid DIFOP calibration from DObject {difop_object!r} under {root}"
+                f"no MSOP packets from DObject {msop_object!r} under {source.label} "
+                f"(log records={len(msop_records)}, "
+                f"host_ticks=[{host_ticks_min}, {host_ticks_max}])"
             )
-        vertical = default_vertical_deg()
-        horizontal = default_horizontal_deg()
-        angle_source = "default_msop_only_vertical_-16_to_16_deg"
-    else:
-        vertical = angles.vertical_deg
-        horizontal = angles.horizontal_deg
-        angle_source = "difop_channel_angles"
 
-    return H32DlogLidarSession(
-        dlog_root=root,
-        msop_object=msop_object,
-        difop_object=difop_object,
-        msop_packets=msop_packets,
-        msop_batch_count=batch_count,
-        difop_record_count=difop_count,
-        angle_source=angle_source,
-        vertical_deg=vertical,
-        horizontal_deg=horizontal,
-        session_id=session_id,
-        lidar_ip=lidar_ip,
-    )
+        if angles is None:
+            if require_difop:
+                raise RuntimeError(
+                    f"no valid DIFOP calibration from DObject {difop_object!r} under {source.label}"
+                )
+            vertical = default_vertical_deg()
+            horizontal = default_horizontal_deg()
+            angle_source = "default_msop_only_vertical_-16_to_16_deg"
+        else:
+            vertical = angles.vertical_deg
+            horizontal = angles.horizontal_deg
+            angle_source = "difop_channel_angles"
+
+        return H32DlogLidarSession(
+            dlog_root=source.label,
+            msop_object=msop_object,
+            difop_object=difop_object,
+            msop_packets=msop_packets,
+            msop_host_utc_ticks=msop_host_utc_ticks,
+            msop_batch_count=batch_count,
+            difop_record_count=difop_count,
+            angle_source=angle_source,
+            vertical_deg=vertical,
+            horizontal_deg=horizontal,
+            session_id=session_id,
+            lidar_ip=lidar_ip,
+            host_ticks_min=host_ticks_min,
+            host_ticks_max=host_ticks_max,
+        )

tools/h32_dlog/timeutil.py

@@ -0,0 +1,56 @@
+"""Wall-clock helpers for Medulla tick filtering.
+
+Two tick conventions appear in this dataset:
+
+- LiDAR DObject ``tic`` / recovered ``indices.log``: ``DateTime.Now.Ticks`` (local)
+- IMU / MSOP payload host receive fields: UTC ``DateTime.UtcNow.Ticks``
+"""
+
+from __future__ import annotations
+
+from datetime import datetime, timedelta, timezone
+
+TICKS_PER_SECOND = 10_000_000
+DOTNET_UNIX_EPOCH_TICKS = 621355968000000000
+
+
+def _parse_local_wall(text: str) -> datetime:
+    normalized = text.strip().replace(" ", "T")
+    if normalized.endswith("Z"):
+        raise ValueError("expected local wall time without Z; got UTC marker")
+    if "+" in normalized[10:]:
+        idx = normalized.find("+", 10)
+        normalized = normalized[:idx]
+    elif normalized.count("-") > 2:
+        # timezone like -08:00 after the date
+        idx = normalized.find("-", 10)
+        if idx > 0 and ":" in normalized[idx + 1 :]:
+            normalized = normalized[:idx]
+    return datetime.fromisoformat(normalized).replace(tzinfo=None)
+
+
+def local_wall_to_dotnet_ticks(text: str) -> int:
+    """Local wall time → ``DateTime.Now.Ticks`` (LiDAR DObject tic)."""
+
+    dt = _parse_local_wall(text)
+    delta = dt - datetime(1, 1, 1)
+    return int(delta.total_seconds() * TICKS_PER_SECOND)
+
+
+def local_wall_to_utc_dotnet_ticks(text: str, *, tz_hours: float = 8.0) -> int:
+    """Local wall time in ``tz_hours`` → UTC ``DateTime.UtcNow.Ticks`` (IMU host)."""
+
+    dt = _parse_local_wall(text).replace(tzinfo=timezone(timedelta(hours=tz_hours)))
+    unix = dt.timestamp()
+    return int(round(unix * TICKS_PER_SECOND)) + DOTNET_UNIX_EPOCH_TICKS
+
+
+def dotnet_ticks_to_local_iso(ticks: int) -> str:
+    dt = datetime(1, 1, 1) + timedelta(microseconds=ticks / 10.0)
+    return dt.isoformat(timespec="milliseconds")
+
+
+def utc_dotnet_ticks_to_unix_s(ticks: int) -> float:
+    """UTC ``DateTime.UtcNow.Ticks`` → Unix seconds."""
+
+    return (float(ticks) - float(DOTNET_UNIX_EPOCH_TICKS)) / float(TICKS_PER_SECOND)

tools/rscap_v2/h32_msop.py

@@ -15,7 +15,7 @@ and horizontal channel offsets default to 0.
 from __future__ import annotations
 
 from dataclasses import dataclass
-from typing import Iterable
+from typing import Iterable, Sequence
 
 import numpy as np
 
@@ -66,6 +66,8 @@ class LidarFrameExport:
     t_start_s: float
     t_end_s: float
     points_xyz: np.ndarray  # (N, 3) metres
+    host_receive_utc_ticks_start: int = 0
+    host_receive_utc_ticks_end: int = 0
 
 
 def decode_packet_points(
@@ -144,6 +146,7 @@ def _block_points(
 def iter_h32_frames_from_packets(
     packets: Iterable[bytes],
     *,
+    host_utc_ticks: Sequence[int] | None = None,
     min_frame_points: int = MIN_FRAME_POINTS_DEFAULT,
     frame_stride: int = 1,
     min_range_m: float = 0.3,
@@ -152,31 +155,49 @@ def iter_h32_frames_from_packets(
     vertical_deg: np.ndarray | None = None,
     horizontal_deg: np.ndarray | None = None,
 ) -> list[LidarFrameExport]:
-    """Assemble raw MSOP packets into frames using the 270°→90° azimuth wrap."""
+    """Assemble raw MSOP packets into frames using the 270°→90° azimuth wrap.
+
+    ``host_utc_ticks`` is optional per-packet ``HostReceiveUtcTicks`` from the
+    MSOP DLog payload (UTC DateTime ticks). When provided, each emitted frame
+    carries host receive start/end ticks from the first/last contributing packet.
+    """
 
     vertical = default_vertical_deg() if vertical_deg is None else np.asarray(vertical_deg, dtype=np.float64)
     horizontal = default_horizontal_deg() if horizontal_deg is None else np.asarray(horizontal_deg, dtype=np.float64)
     if vertical.shape != (CHANNELS,) or horizontal.shape != (CHANNELS,):
         raise ValueError(f"vertical/horizontal must have shape ({CHANNELS},)")
 
+    packet_list = list(packets)
+    host_list = list(host_utc_ticks) if host_utc_ticks is not None else None
+    if host_list is not None and len(host_list) != len(packet_list):
+        raise ValueError(
+            f"host_utc_ticks length {len(host_list)} != packets length {len(packet_list)}"
+        )
+
     frames: list[LidarFrameExport] = []
     point_chunks: list[np.ndarray] = []
     t_start: float | None = None
     t_end: float | None = None
+    host_start: int | None = None
+    host_end: int | None = None
     prev_az: float | None = None
     kept = 0
     stride = max(1, int(frame_stride))
 
     def emit() -> None:
-        nonlocal point_chunks, t_start, t_end, kept
+        nonlocal point_chunks, t_start, t_end, host_start, host_end, kept
         if not point_chunks or t_start is None or t_end is None:
             point_chunks = []
             t_start = t_end = None
+            host_start = host_end = None
             return
         points = np.vstack(point_chunks)
         point_chunks = []
         start_s, end_s = t_start, t_end
+        h0 = int(host_start or 0)
+        h1 = int(host_end or 0)
         t_start = t_end = None
+        host_start = host_end = None
         if points.shape[0] < min_frame_points:
             return
         if kept % stride != 0:
@@ -188,12 +209,21 @@ def iter_h32_frames_from_packets(
             points = points[select]
         if end_s <= start_s:
             end_s = start_s + 0.1
-        frames.append(LidarFrameExport(t_start_s=start_s, t_end_s=end_s, points_xyz=points))
+        frames.append(
+            LidarFrameExport(
+                t_start_s=start_s,
+                t_end_s=end_s,
+                points_xyz=points,
+                host_receive_utc_ticks_start=h0,
+                host_receive_utc_ticks_end=h1,
+            )
+        )
 
-    for packet in packets:
+    for index, packet in enumerate(packet_list):
         if len(packet) != PACKET_LENGTH:
             continue
         packet_t = device_timestamp_ms(packet) * 1e-3
+        packet_host = int(host_list[index]) if host_list is not None else 0
         unit = distance_unit_mm(packet)
         idx = DATA_START
         for _block in range(BLOCKS):
@@ -216,7 +246,9 @@ def iter_h32_frames_from_packets(
             if pts.shape[0]:
                 if t_start is None:
                     t_start = packet_t
+                    host_start = packet_host
                 t_end = packet_t
+                host_end = packet_host
                 point_chunks.append(pts)
             idx += BLOCK_LENGTH
 
@@ -237,8 +269,11 @@ def iter_h32_frames(
 ) -> list[LidarFrameExport]:
     """Assemble MSOP packets from a V2 .rscap capture into frames."""
 
+    packets = [chunk.raw for chunk in capture.chunks]
+    host_ticks = [chunk.receive_utc_ticks for chunk in capture.chunks]
     return iter_h32_frames_from_packets(
-        (chunk.raw for chunk in capture.chunks),
+        packets,
+        host_utc_ticks=host_ticks,
         min_frame_points=min_frame_points,
         frame_stride=frame_stride,
         min_range_m=min_range_m,

tools/rscap_v2/hi13_imu.py

@@ -0,0 +1,136 @@
+"""Decode Hipnuc / HI13 (HI91/HI92) IMU frames from a V2 .rscap capture.
+
+Matches ``EcarSensorMinimal/RawSerialImu/Hi13Protocol.cs``:
+sync ``5A A5``, CRC16 over header[0:4]+payload, tag ``0x91`` / ``0x92``.
+
+HI91 (preferred for calibration):
+- accel: float32 in g → m/s² (* 9.80665)
+- gyro: float32 in deg/s → rad/s
+- device time: uint32 ms at frame offset 14 → ``t_s = ms * 1e-3``
+"""
+
+from __future__ import annotations
+
+import struct
+
+import numpy as np
+
+from .capture_format_v2 import CaptureFile
+from .n300_imu import ImuSample, samples_to_arrays
+
+G0 = 9.80665
+DEG2RAD = np.pi / 180.0
+
+
+def crc16_hi13(frame: bytes, payload_length: int) -> int:
+    crc = 0
+    for value in frame[:4]:
+        crc = _update_crc16(crc, value)
+    for value in frame[6 : 6 + payload_length]:
+        crc = _update_crc16(crc, value)
+    return crc & 0xFFFF
+
+
+def _update_crc16(crc: int, value: int) -> int:
+    crc ^= (value & 0xFF) << 8
+    for _ in range(8):
+        if crc & 0x8000:
+            crc = ((crc << 1) ^ 0x1021) & 0xFFFF
+        else:
+            crc = (crc << 1) & 0xFFFF
+    return crc
+
+
+def parse_hi91_frame(raw: bytes) -> tuple[tuple[float, float, float], tuple[float, float, float], int] | None:
+    """Return (gyro_rad_s, accel_m_s2, device_timestamp_ms) for a CRC-valid HI91 frame."""
+
+    if len(raw) < 6 + 76:
+        return None
+    payload_length = raw[2] | (raw[3] << 8)
+    if payload_length < 76 or len(raw) < 6 + payload_length:
+        return None
+    if raw[6] != 0x91:
+        return None
+    expected = raw[4] | (raw[5] << 8)
+    if crc16_hi13(raw, payload_length) != expected:
+        return None
+    device_ms = struct.unpack_from("<I", raw, 14)[0]
+    ax, ay, az = struct.unpack_from("<fff", raw, 18)
+    gx, gy, gz = struct.unpack_from("<fff", raw, 30)
+    gyro = (gx * DEG2RAD, gy * DEG2RAD, gz * DEG2RAD)
+    accel = (ax * G0, ay * G0, az * G0)
+    return gyro, accel, int(device_ms)
+
+
+def iter_hi13_imu_samples(
+    capture: CaptureFile,
+    *,
+    host_utc_ticks_min: int | None = None,
+    host_utc_ticks_max: int | None = None,
+) -> list[ImuSample]:
+    """Return CRC-valid HI91 samples sorted by device timestamp.
+
+    Streams chunk-by-chunk (no giant join) and can skip whole chunks outside the
+    host UTC receive window before parsing.
+    """
+
+    samples: list[ImuSample] = []
+    carry = b""
+    for chunk in capture.chunks:
+        if host_utc_ticks_min is not None and chunk.receive_utc_ticks < host_utc_ticks_min:
+            carry = b""
+            continue
+        if host_utc_ticks_max is not None and chunk.receive_utc_ticks > host_utc_ticks_max:
+            # chunks are time-ordered; remaining ones are later
+            if chunk.receive_utc_ticks > host_utc_ticks_max:
+                break
+        stream = carry + chunk.raw
+        cursor = 0
+        while cursor + 6 < len(stream):
+            sync = stream.find(b"\x5A\xA5", cursor)
+            if sync < 0:
+                carry = b""
+                break
+            if sync + 6 > len(stream):
+                carry = stream[sync:]
+                break
+            payload_length = stream[sync + 2] | (stream[sync + 3] << 8)
+            if payload_length < 1 or payload_length > 512:
+                cursor = sync + 1
+                continue
+            end = sync + 6 + payload_length
+            if end > len(stream):
+                carry = stream[sync:]
+                break
+            parsed = parse_hi91_frame(stream[sync:end])
+            cursor = end
+            if parsed is None:
+                continue
+            gyro, accel, device_ms = parsed
+            host_ticks = chunk.receive_utc_ticks
+            if host_utc_ticks_min is not None and host_ticks < host_utc_ticks_min:
+                continue
+            if host_utc_ticks_max is not None and host_ticks > host_utc_ticks_max:
+                continue
+            samples.append(
+                ImuSample(
+                    t_s=float(device_ms) * 1e-3,
+                    gyro_rad_s=gyro,
+                    accel_m_s2=accel,
+                    host_receive_utc_ticks=host_ticks,
+                    device_timestamp_us=int(device_ms) * 1000,
+                )
+            )
+        else:
+            carry = b""
+    samples.sort(key=lambda sample: (sample.t_s, sample.device_timestamp_us))
+    return samples
+
+
+__all__ = [
+    "ImuSample",
+    "crc16_hi13",
+    "iter_hi13_imu_samples",
+    "parse_hi91_frame",
+    "samples_to_arrays",
+]

tools/run_joint_three_windows_full_se3.py

@@ -0,0 +1,50 @@
+#!/usr/bin/env python3
+"""Joint full_se3 on the three host-aligned priority windows."""
+
+from __future__ import annotations
+
+import sys
+from pathlib import Path
+
+ROOT = Path(__file__).resolve().parents[1]
+if str(ROOT) not in sys.path:
+    sys.path.insert(0, str(ROOT))
+
+from imu_lidar.cli import main as cli_main
+
+ALIGNED = Path(r"D:\data\calibration_usable_20260808\sessions_v1_host_aligned")
+SESSIONS = [
+    "priority_174005_174515",
+    "priority_174905_175450",
+    "priority_175910_180530",
+]
+
+
+def main() -> int:
+    out = ALIGNED / "joint_full_se3"
+    argv = [
+        "run",
+        "--vehicle-config",
+        str(ROOT / "config" / "vehicle_hi13_h32_20260808.yaml"),
+        "--output",
+        str(out),
+        "--mode",
+        "full_se3",
+        "--time-offset-search-s",
+        "0.5",
+        "--min-pair-rotation-deg",
+        "2.0",
+    ]
+    for name in SESSIONS:
+        session = ALIGNED / name
+        if not (session / "imu.csv").is_file() or not (session / "lidar" / "frames_index.csv").is_file():
+            raise SystemExit(f"missing host-aligned session: {session}")
+        argv.extend(["--session-id", name])
+        argv.extend(["--imu", str(session / "imu.csv")])
+        argv.extend(["--lidar", str(session / "lidar")])
+    print("argv:", " ".join(argv), flush=True)
+    return cli_main(argv)
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

tools/run_priority_windows_calibration.py

@@ -0,0 +1,292 @@
+#!/usr/bin/env python3
+"""Align HI13/H32 via host-UTC bridge, then run rotation_only.
+
+Device clocks (HI13 boot ms vs H32 absolute) must NOT be forced to share a
+first-sample epoch. Instead map each LiDAR frame onto the IMU device timeline
+by interpolating IMU device time at the frame's MSOP HostReceiveUtcTicks.
+Optional |ω| correlation then refines residual host/path delay.
+"""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import json
+import shutil
+import sys
+from pathlib import Path
+
+import numpy as np
+
+ROOT = Path(__file__).resolve().parents[1]
+if str(ROOT) not in sys.path:
+    sys.path.insert(0, str(ROOT))
+
+from imu_lidar.cli import main as cli_main
+from imu_lidar.geometry import so3_log
+from imu_lidar.imu_io import load_imu_samples
+from imu_lidar.lidar_io import load_lidar_frames
+from imu_lidar.registration import estimate_frame_rotations
+from imu_lidar.time_offset import _correlate_offset, _magnitude_series
+
+SESSIONS = [
+    "priority_174005_174515",
+    "priority_174905_175450",
+    "priority_175910_180530",
+]
+
+
+def _read_imu_host_table(imu_csv: Path) -> tuple[np.ndarray, np.ndarray]:
+    rows = list(csv.DictReader(imu_csv.open(encoding="utf-8")))
+    if not rows:
+        raise RuntimeError(f"empty IMU csv: {imu_csv}")
+    if "t_host_utc_s" not in rows[0] or not rows[0].get("t_host_utc_s"):
+        raise RuntimeError(
+            f"{imu_csv} missing t_host_utc_s; re-export with HostReceiveUtcTicks support"
+        )
+    t_dev = np.asarray([float(row["t"]) for row in rows], dtype=np.float64)
+    t_host = np.asarray([float(row["t_host_utc_s"]) for row in rows], dtype=np.float64)
+    order = np.argsort(t_host)
+    return t_host[order], t_dev[order]
+
+
+def _imu_device_at_host(t_host_query: np.ndarray, imu_host: np.ndarray, imu_dev: np.ndarray) -> np.ndarray:
+    """Map host UTC seconds → IMU device seconds (linear interp, edge clamp)."""
+
+    return np.interp(t_host_query, imu_host, imu_dev)
+
+
+def rewrite_lidar_index_host_bridge(
+    src_index: Path,
+    dst_index: Path,
+    *,
+    imu_host: np.ndarray,
+    imu_dev: np.ndarray,
+    residual_delta_s: float = 0.0,
+) -> dict:
+    """Rewrite LiDAR times onto IMU device clock via host UTC bridge.
+
+    For each frame:
+      t_host_mid = mid of MSOP host receive window
+      t_imu_mid  = interp(IMU device @ t_host_mid) + residual_delta
+      keep device duration: t_start/t_end centered on t_imu_mid
+    """
+
+    rows = list(csv.DictReader(src_index.open(encoding="utf-8")))
+    if not rows:
+        raise RuntimeError(f"empty frames_index: {src_index}")
+    if "t_host_utc_s" not in rows[0]:
+        raise RuntimeError(
+            f"{src_index} missing t_host_utc_s; re-export DLog with MSOP HostReceiveUtcTicks"
+        )
+
+    dst_index.parent.mkdir(parents=True, exist_ok=True)
+    offsets: list[float] = []
+    with dst_index.open("w", newline="", encoding="utf-8") as handle:
+        writer = csv.writer(handle)
+        writer.writerow(["frame_id", "filename", "t_start", "t_end"])
+        for row in rows:
+            t0 = float(row["t_start"])
+            t1 = float(row["t_end"])
+            host0 = row.get("t_host_utc_s") or ""
+            host1 = row.get("t_host_utc_end_s") or ""
+            if not host0:
+                raise RuntimeError(f"frame {row.get('frame_id')} missing t_host_utc_s")
+            h0 = float(host0)
+            h1 = float(host1) if host1 else h0
+            host_mid = 0.5 * (h0 + h1)
+            imu_mid = float(_imu_device_at_host(np.asarray([host_mid]), imu_host, imu_dev)[0])
+            imu_mid += residual_delta_s
+            duration = max(t1 - t0, 1e-3)
+            new0 = imu_mid - 0.5 * duration
+            new1 = imu_mid + 0.5 * duration
+            offsets.append(imu_mid - 0.5 * (t0 + t1))
+            writer.writerow(
+                [
+                    row["frame_id"],
+                    row["filename"],
+                    f"{new0:.9f}",
+                    f"{new1:.9f}",
+                ]
+            )
+    arr = np.asarray(offsets, dtype=np.float64)
+    return {
+        "frames": len(offsets),
+        "bridge_offset_median_s": float(np.median(arr)),
+        "bridge_offset_mean_s": float(np.mean(arr)),
+        "bridge_offset_std_s": float(np.std(arr)),
+        "bridge_offset_min_s": float(np.min(arr)),
+        "bridge_offset_max_s": float(np.max(arr)),
+        "residual_delta_s": float(residual_delta_s),
+    }
+
+
+def estimate_residual_delta(session_dir: Path, *, search_s: float = 5.0) -> tuple[float, float]:
+    imu = load_imu_samples(session_dir / "imu.csv")
+    frames = load_lidar_frames(session_dir / "lidar")
+    # Short pairs only — large stride anti-correlates with IMU |gyro|.
+    stride = 1 if len(frames) < 80 else 2
+    rotations, pair_times = estimate_frame_rotations(frames, stride=stride)
+    if len(rotations) < 8:
+        rotations, pair_times = estimate_frame_rotations(frames, stride=1)
+    lidar_t = []
+    lidar_w = []
+    for (t_a, t_b), rotation in zip(pair_times, rotations):
+        dt_pair = max(t_b - t_a, 1e-3)
+        omega = so3_log(rotation) / dt_pair
+        lidar_t.append(0.5 * (t_a + t_b))
+        lidar_w.append(omega)
+    imu_t, imu_mag = _magnitude_series(imu.t_s, imu.gyro_rad_s)
+    lidar_t_arr, lidar_mag = _magnitude_series(np.asarray(lidar_t), np.asarray(lidar_w))
+    delta, peak = _correlate_offset(
+        imu_t,
+        imu_mag,
+        lidar_t_arr,
+        lidar_mag,
+        search_s=search_s,
+        sample_hz=20.0,
+    )
+    return float(delta), float(peak)
+
+
+def align_session(src: Path, dst: Path, *, residual_search_s: float = 5.0) -> dict:
+    if dst.exists():
+        shutil.rmtree(dst)
+    dst.mkdir(parents=True)
+    shutil.copy2(src / "imu.csv", dst / "imu.csv")
+    shutil.copytree(src / "lidar" / "frames", dst / "lidar" / "frames")
+
+    imu_host, imu_dev = _read_imu_host_table(src / "imu.csv")
+    bridge = rewrite_lidar_index_host_bridge(
+        src / "lidar" / "frames_index.csv",
+        dst / "lidar" / "frames_index.csv",
+        imu_host=imu_host,
+        imu_dev=imu_dev,
+        residual_delta_s=0.0,
+    )
+
+    residual_delta, residual_peak = estimate_residual_delta(dst, search_s=residual_search_s)
+    # Only apply residual when correlation is clearly positive; otherwise the
+    # host-UTC bridge alone is the trusted alignment (weak peaks are noise).
+    apply_residual = residual_peak >= 0.5 and abs(residual_delta) <= residual_search_s
+    applied = float(residual_delta) if apply_residual else 0.0
+    if apply_residual:
+        bridge = rewrite_lidar_index_host_bridge(
+            src / "lidar" / "frames_index.csv",
+            dst / "lidar" / "frames_index.csv",
+            imu_host=imu_host,
+            imu_dev=imu_dev,
+            residual_delta_s=applied,
+        )
+
+    meta = {
+        "source": str(src),
+        "aligned": str(dst),
+        "method": "host_utc_bridge",
+        "imu_host_span_s": [float(imu_host[0]), float(imu_host[-1])],
+        "imu_device_span_s": [float(imu_dev[0]), float(imu_dev[-1])],
+        "bridge": bridge,
+        "residual_delta_s": residual_delta,
+        "residual_peak": residual_peak,
+        "residual_applied_s": applied,
+        "residual_applied": apply_residual,
+        "note": (
+            "LiDAR t_* rewritten onto IMU device clock via MSOP/IMU HostReceiveUtc; "
+            "not first-device-sample coincidence. residual |omega| shift applied only if peak>=0.5."
+        ),
+    }
+    (dst / "align_meta.json").write_text(
+        json.dumps(meta, indent=2, ensure_ascii=False) + "\n", encoding="utf-8"
+    )
+    return meta
+
+
+def run_one(session_dir: Path, vehicle: Path, search_s: float) -> dict:
+    output = session_dir / "out"
+    if output.exists():
+        shutil.rmtree(output)
+    argv = [
+        "run",
+        "--session-id",
+        session_dir.name,
+        "--imu",
+        str(session_dir / "imu.csv"),
+        "--lidar",
+        str(session_dir / "lidar"),
+        "--vehicle-config",
+        str(vehicle),
+        "--output",
+        str(output),
+        "--mode",
+        "rotation_only",
+        "--time-offset-search-s",
+        str(search_s),
+        "--min-pair-rotation-deg",
+        "2.0",
+    ]
+    code = cli_main(argv)
+    summary_path = output / "summary.json"
+    summary = {}
+    if summary_path.is_file():
+        summary = json.loads(summary_path.read_text(encoding="utf-8"))
+    t_block = summary.get("T_IMU_lidar") or {}
+    return {
+        "session": session_dir.name,
+        "exit_code": code,
+        "status": summary.get("status"),
+        "message": summary.get("message"),
+        "time_offset_s": summary.get("time_offset_s"),
+        "rotation_deg": t_block.get("rotation_deg") if isinstance(t_block, dict) else None,
+        "summary": str(summary_path) if summary_path.is_file() else None,
+    }
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument(
+        "--sessions-root",
+        type=Path,
+        default=Path(r"D:\data\calibration_usable_20260808\sessions_v1"),
+    )
+    parser.add_argument(
+        "--aligned-root",
+        type=Path,
+        default=Path(r"D:\data\calibration_usable_20260808\sessions_v1_host_aligned"),
+    )
+    parser.add_argument(
+        "--vehicle-config",
+        type=Path,
+        default=ROOT / "config" / "vehicle_hi13_h32_20260808.yaml",
+    )
+    parser.add_argument(
+        "--residual-search-s",
+        type=float,
+        default=5.0,
+        help="|ω| residual search after host bridge (seconds)",
+    )
+    parser.add_argument("--time-offset-search-s", type=float, default=1.0)
+    args = parser.parse_args()
+
+    results = []
+    for name in SESSIONS:
+        src = args.sessions_root / name
+        if not src.is_dir():
+            raise SystemExit(f"missing session: {src}")
+        aligned = args.aligned_root / name
+        print(f"=== align {name} ===", flush=True)
+        meta = align_session(src, aligned, residual_search_s=args.residual_search_s)
+        print(json.dumps(meta, ensure_ascii=False, indent=2), flush=True)
+        print(f"=== calibrate {name} ===", flush=True)
+        result = run_one(aligned, args.vehicle_config, args.time_offset_search_s)
+        results.append({"align": meta, **result})
+        print(json.dumps(result, ensure_ascii=False, indent=2), flush=True)
+
+    manifest = args.aligned_root / "calibration_manifest.json"
+    args.aligned_root.mkdir(parents=True, exist_ok=True)
+    manifest.write_text(json.dumps(results, ensure_ascii=False, indent=2) + "\n", encoding="utf-8")
+    print(f"manifest: {manifest}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

tools/visualize_pair_3d.py

@@ -9,6 +9,8 @@ Modes (keyboard):
   2          IMU prediction with X=I  (B_pred = A)
   3          LiDAR registration B     (reference)
   4          calibrated prediction    B_pred = X^{-1} A X
+             (rotation_only runs default to R conjug + t_B so bad IMU Δp
+              does not dominate the overlay)
   N / ]      next motion pair
   P / [      previous motion pair
   Q / Esc    exit
@@ -20,10 +22,16 @@ from __future__ import annotations
 
 import argparse
 import json
+import sys
 from pathlib import Path
+from typing import Any
 
 import numpy as np
 
+ROOT = Path(__file__).resolve().parents[1]
+if str(ROOT) not in sys.path:
+    sys.path.insert(0, str(ROOT))
+
 from imu_lidar.geometry import (
     inverse_transform,
     make_transform,
@@ -34,12 +42,39 @@ from imu_lidar.geometry import (
 from imu_lidar.imu_io import load_imu_samples
 from imu_lidar.imu_preintegration import preintegrate_imu
 from imu_lidar.keyframes import build_keyframes
-from imu_lidar.lidar_io import load_lidar_frames
+from imu_lidar.lidar_io import load_lidar_frame_at, load_lidar_frames
 from imu_lidar.motion_pairs import build_motion_pairs
+from imu_lidar.motion_pairs_io import (
+    load_motion_pairs,
+    pairs_for_session,
+    resolve_motion_pairs_path,
+)
 from imu_lidar.registration import register_lidar_pair
 from imu_lidar.time_offset import lidar_time_to_imu_time
 
 
+class _LazyFrameStore:
+    """Load NPZ frames on demand; indices match mid-time-sorted ``load_lidar_frames``."""
+
+    def __init__(self, lidar_dir: Path, *, max_cached: int = 16):
+        self.lidar_dir = Path(lidar_dir)
+        self.max_cached = max_cached
+        self._cache: dict[int, object] = {}
+        self._order: list[int] = []
+
+    def __getitem__(self, index: int):
+        index = int(index)
+        if index in self._cache:
+            return self._cache[index]
+        frame = load_lidar_frame_at(self.lidar_dir, index)
+        self._cache[index] = frame
+        self._order.append(index)
+        while len(self._order) > self.max_cached:
+            old = self._order.pop(0)
+            self._cache.pop(old, None)
+        return frame
+
+
 COLORS = {
     "target": [0.10, 0.65, 1.00],
     "source": [1.00, 0.35, 0.05],
@@ -53,24 +88,45 @@ MODE_NAMES = (
 )
 
 
-def _load_extrinsic(summary_path: Path) -> tuple[np.ndarray, float, np.ndarray]:
+def _load_extrinsic(summary_path: Path) -> tuple[np.ndarray, float, np.ndarray, dict[str, Any]]:
     summary = json.loads(summary_path.read_text(encoding="utf-8"))
     t_block = summary.get("T_IMU_lidar")
+    meta: dict[str, Any] = {
+        "rotation_only": False,
+        "translation_accepted": False,
+        "status": str(summary.get("status") or ""),
+    }
     if t_block is None:
         matrix = summary.get("matrix")
         if matrix is not None:
-            return np.asarray(matrix, dtype=float), 0.0, np.zeros(3)
+            return np.asarray(matrix, dtype=float), 0.0, np.zeros(3), meta
         raise ValueError(f"no T_IMU_lidar in {summary_path}")
     t_mat = np.asarray(t_block["matrix"], dtype=float)
     dt = float(summary.get("time_offset_s") or 0.0)
     session = (summary.get("details") or {}).get("sessions", [{}])[0]
+    joint = session.get("joint") or {}
     bias = np.asarray(
         (session.get("imu_audit") or {}).get("gyro_bias_rad_s")
-        or (session.get("joint") or {}).get("gyro_bias_rad_s")
+        or joint.get("gyro_bias_rad_s")
         or [0.0, 0.0, 0.0],
         dtype=float,
     ).reshape(3)
-    return t_mat, dt, bias
+    status = str(summary.get("status") or "")
+    translation_accepted = bool(
+        joint.get("translation_accepted")
+        or (summary.get("details") or {}).get("joint", {}).get("translation_accepted")
+    )
+    rotation_only = ("rotation_only" in status) or (
+        not translation_accepted and float(np.linalg.norm(t_mat[:3, 3])) < 1e-9
+    )
+    meta.update(
+        {
+            "rotation_only": rotation_only,
+            "translation_accepted": translation_accepted,
+            "status": status,
+        }
+    )
+    return t_mat, dt, bias, meta
 
 
 def _delta_components(reference: np.ndarray, candidate: np.ndarray) -> dict:
@@ -169,16 +225,42 @@ def _pair_from_indices(
     return frame_i, frame_j, a, reg.transform
 
 
-def _transforms_for_pair(x: np.ndarray, a_ij: np.ndarray, b_gicp: np.ndarray) -> dict[str, np.ndarray]:
+def _transforms_for_pair(
+    x: np.ndarray,
+    a_ij: np.ndarray,
+    b_gicp: np.ndarray,
+    *,
+    mode4_translation: str = "imu",
+) -> dict[str, np.ndarray]:
+    """Build overlay transforms.
+
+    ``mode4_translation``:
+      - ``imu``: full SE3 conjug ``X^{-1} A X`` (needs trustworthy IMU Δp)
+      - ``gicp``: rotation conjug only; translation taken from LiDAR B
+        (correct check for rotation_only calibrations)
+    """
+
+    calibrated = inverse_transform(x) @ a_ij @ x
+    if mode4_translation == "gicp":
+        calibrated = make_transform(b_gicp[:3, 3], calibrated[:3, :3])
+    elif mode4_translation != "imu":
+        raise ValueError(f"unknown mode4_translation={mode4_translation!r}")
     return {
         MODE_NAMES[0]: np.eye(4),
         MODE_NAMES[1]: a_ij.copy(),
         MODE_NAMES[2]: b_gicp.copy(),
-        MODE_NAMES[3]: inverse_transform(x) @ a_ij @ x,
+        MODE_NAMES[3]: calibrated,
     }
 
 
-def _resolve_pair(frames, pairs, pair_index: int, x: np.ndarray):
+def _resolve_pair(
+    frames,
+    pairs,
+    pair_index: int,
+    x: np.ndarray,
+    *,
+    mode4_translation: str = "imu",
+):
     pair = pairs[pair_index]
     frame_i = frames[pair.i]
     frame_j = frames[pair.j]
@@ -190,7 +272,9 @@ def _resolve_pair(frames, pairs, pair_index: int, x: np.ndarray):
         pair.t_B_m if pair.t_B_m is not None else np.zeros(3),
         pair.R_B,
     )
-    transforms = _transforms_for_pair(x, a_ij, b_gicp)
+    transforms = _transforms_for_pair(
+        x, a_ij, b_gicp, mode4_translation=mode4_translation
+    )
     label = (
         f"pair {pair_index + 1}/{len(pairs)}  "
         f"frames {pair.i} <- {pair.j}  "
@@ -255,12 +339,15 @@ def _run_gui(
     start_index: int,
     voxel: float,
     fixed_single_pair: tuple | None,
+    mode4_translation: str = "imu",
 ) -> None:
     import open3d as o3d
 
     if fixed_single_pair is not None:
         frame_i, frame_j, a_ij, b_gicp = fixed_single_pair
-        transforms = _transforms_for_pair(x, a_ij, b_gicp)
+        transforms = _transforms_for_pair(
+            x, a_ij, b_gicp, mode4_translation=mode4_translation
+        )
         label = f"fixed frames (no pair switching)"
         pair_index = 0
         n_pairs = 1
@@ -270,7 +357,7 @@ def _run_gui(
         n_pairs = len(pairs)
         use_list = True
         frame_i, frame_j, a_ij, b_gicp, transforms, label = _resolve_pair(
-            frames, pairs, pair_index, x
+            frames, pairs, pair_index, x, mode4_translation=mode4_translation
         )
 
     viewer = o3d.visualization.VisualizerWithKeyCallback()
@@ -311,7 +398,7 @@ def _run_gui(
             return
         new_index = int(new_index) % n_pairs
         frame_i, frame_j, _a, b_gicp, transforms, label = _resolve_pair(
-            frames, pairs, new_index, x
+            frames, pairs, new_index, x, mode4_translation=mode4_translation
         )
         state["pair_index"] = new_index
         state["transforms"] = transforms
@@ -357,13 +444,34 @@ def _run_gui(
 def main(argv: list[str] | None = None) -> int:
     parser = argparse.ArgumentParser(description=__doc__)
     parser.add_argument("--lidar", required=True, type=Path, help="LiDAR session directory")
-    parser.add_argument("--imu", required=True, type=Path, help="IMU CSV")
+    parser.add_argument(
+        "--imu",
+        type=Path,
+        default=None,
+        help="IMU CSV (only needed when rebuilding pairs without motion_pairs.json)",
+    )
     parser.add_argument(
         "--summary",
         required=True,
         type=Path,
         help="summary.json (or T_IMU_lidar.json) from a calibration run",
     )
+    parser.add_argument(
+        "--motion-pairs",
+        type=Path,
+        default=None,
+        help="Cached motion_pairs.json (default: next to --summary)",
+    )
+    parser.add_argument(
+        "--session-id",
+        default=None,
+        help="Session id inside multi-session motion_pairs.json",
+    )
+    parser.add_argument(
+        "--rebuild-pairs",
+        action="store_true",
+        help="Ignore cache and rebuild pairs from IMU/LiDAR (slow)",
+    )
     parser.add_argument("--pair-index", type=int, default=0, help="Starting motion-pair index")
     parser.add_argument("--frame-i", type=int, default=None, help="Optional explicit frame index i")
     parser.add_argument("--frame-j", type=int, default=None, help="Optional explicit frame index j")
@@ -381,45 +489,104 @@ def main(argv: list[str] | None = None) -> int:
         action="store_true",
         help="Skip Open3D window (use with --save-png)",
     )
+    parser.add_argument(
+        "--mode4-translation",
+        choices=("auto", "gicp", "imu"),
+        default="auto",
+        help=(
+            "Mode-4 translation source: gicp=R conjug + t_B (rotation check); "
+            "imu=full X^-1 A X; auto=gicp for rotation_only summaries"
+        ),
+    )
     args = parser.parse_args(argv)
 
-    x, delta_t_s, gyro_bias = _load_extrinsic(args.summary)
-    frames, imu, keyframes, pair_set = _build_pair_list(
-        lidar_dir=args.lidar,
-        imu_path=args.imu,
-        delta_t_s=delta_t_s,
-        gyro_bias=gyro_bias,
-        min_rotation_deg=args.min_pair_rotation_deg,
-        min_translation_m=args.min_pair_translation_m,
-    )
+    x, delta_t_s, gyro_bias, extr_meta = _load_extrinsic(args.summary)
+    if args.mode4_translation == "auto":
+        mode4_translation = "gicp" if extr_meta.get("rotation_only") else "imu"
+    else:
+        mode4_translation = args.mode4_translation
+    cache_path = args.motion_pairs or resolve_motion_pairs_path(args.summary)
+    use_cache = (not args.rebuild_pairs) and cache_path is not None and args.frame_i is None
 
+    frames = None
+    pairs: tuple = ()
     fixed_single_pair = None
-    if args.frame_i is not None and args.frame_j is not None:
-        frame_i, frame_j, a_ij, b_gicp = _pair_from_indices(
-            frames,
-            imu,
-            i=args.frame_i,
-            j=args.frame_j,
+    label = ""
+    b_gicp = np.eye(4)
+    transforms: dict[str, np.ndarray] = {}
+    frame_i = frame_j = None
+
+    if use_cache:
+        payload = load_motion_pairs(cache_path)
+        pair_list = pairs_for_session(payload, args.session_id)
+        if not pair_list:
+            raise SystemExit(f"no pairs in cache: {cache_path}")
+        if not 0 <= args.pair_index < len(pair_list):
+            raise SystemExit(
+                f"pair-index {args.pair_index} outside [0, {len(pair_list) - 1}] "
+                f"({len(pair_list)} pairs in cache)"
+            )
+        frames = _LazyFrameStore(args.lidar)
+        pairs = tuple(pair_list)
+        frame_i, frame_j, a_ij, b_gicp, transforms, label = _resolve_pair(
+            frames, pairs, args.pair_index, x, mode4_translation=mode4_translation
+        )
+        print(f"loaded {len(pairs)} cached pairs from {cache_path}")
+    else:
+        if args.imu is None:
+            raise SystemExit("--imu is required when motion_pairs.json is missing (or use --rebuild-pairs with --imu)")
+        frames, imu, keyframes, pair_set = _build_pair_list(
+            lidar_dir=args.lidar,
+            imu_path=args.imu,
             delta_t_s=delta_t_s,
             gyro_bias=gyro_bias,
+            min_rotation_deg=args.min_pair_rotation_deg,
+            min_translation_m=args.min_pair_translation_m,
         )
-        transforms = _transforms_for_pair(x, a_ij, b_gicp)
-        label = f"frames {args.frame_i} <- {args.frame_j}"
-        fixed_single_pair = (frame_i, frame_j, a_ij, b_gicp)
-        pairs = ()
-    else:
-        if not pair_set.pairs:
-            raise SystemExit("no motion pairs rebuilt; loosen min-pair thresholds or check data")
-        if not 0 <= args.pair_index < len(pair_set.pairs):
-            raise SystemExit(
-                f"pair-index {args.pair_index} outside [0, {len(pair_set.pairs) - 1}] "
-                f"({len(pair_set.pairs)} pairs available)"
+        if args.frame_i is not None and args.frame_j is not None:
+            frame_i, frame_j, a_ij, b_gicp = _pair_from_indices(
+                frames,
+                imu,
+                i=args.frame_i,
+                j=args.frame_j,
+                delta_t_s=delta_t_s,
+                gyro_bias=gyro_bias,
             )
-        pairs = pair_set.pairs
-        frame_i, frame_j, a_ij, b_gicp, transforms, label = _resolve_pair(
-            frames, pairs, args.pair_index, x
+            transforms = _transforms_for_pair(
+                x, a_ij, b_gicp, mode4_translation=mode4_translation
+            )
+            label = f"frames {args.frame_i} <- {args.frame_j}"
+            fixed_single_pair = (frame_i, frame_j, a_ij, b_gicp)
+            pairs = ()
+        else:
+            if not pair_set.pairs:
+                raise SystemExit("no motion pairs rebuilt; loosen min-pair thresholds or check data")
+            if not 0 <= args.pair_index < len(pair_set.pairs):
+                raise SystemExit(
+                    f"pair-index {args.pair_index} outside [0, {len(pair_set.pairs) - 1}] "
+                    f"({len(pair_set.pairs)} pairs available)"
+                )
+            pairs = pair_set.pairs
+            frame_i, frame_j, a_ij, b_gicp, transforms, label = _resolve_pair(
+                frames, pairs, args.pair_index, x, mode4_translation=mode4_translation
+            )
+            print(f"rebuilt {len(pairs)} pairs from {len(keyframes.indices)} keyframes")
+
+    print(
+        f"mode4 translation={mode4_translation}  "
+        f"(status={extr_meta.get('status') or 'n/a'}, "
+        f"rotation_only={bool(extr_meta.get('rotation_only'))})"
+    )
+    if mode4_translation == "gicp":
+        print(
+            "note: mode4 uses R conjug + t_B; IMU Δp is ignored "
+            "(typical for rotation_only — raw Δp often has large Z drift)."
+        )
+    if mode4_translation == "imu":
+        print(
+            "note: mode4 uses full X^-1 A X. If clouds stack vertically, "
+            "IMU Δp is likely bad; retry with --mode4-translation gicp."
         )
-        print(f"rebuilt {len(pairs)} pairs from {len(keyframes.indices)} keyframes")
 
     if args.save_png is not None:
         _print_pair_header(label, b_gicp, transforms)
@@ -444,6 +611,7 @@ def main(argv: list[str] | None = None) -> int:
         start_index=args.pair_index,
         voxel=args.voxel,
         fixed_single_pair=fixed_single_pair,
+        mode4_translation=mode4_translation,
     )
     return 0