Get new device to head transform based on good segments. Segments starting with "BAD" annotations are not included for calculating the mean head position. Parameters ---------- raw : instance of Raw | list of Raw Data to compute head position. Can be a list containing m
(raw, pos, *, verbose=None)
| 301 | |
| 302 | @verbose |
| 303 | def compute_average_dev_head_t(raw, pos, *, verbose=None): |
| 304 | """Get new device to head transform based on good segments. |
| 305 | |
| 306 | Segments starting with "BAD" annotations are not included for calculating |
| 307 | the mean head position. |
| 308 | |
| 309 | Parameters |
| 310 | ---------- |
| 311 | raw : instance of Raw | list of Raw |
| 312 | Data to compute head position. Can be a list containing multiple raw |
| 313 | instances. |
| 314 | pos : array, shape (N, 10) | list of ndarray |
| 315 | The position and quaternion parameters from cHPI fitting. Can be |
| 316 | a list containing multiple position arrays, one per raw instance passed. |
| 317 | %(verbose)s |
| 318 | |
| 319 | Returns |
| 320 | ------- |
| 321 | dev_head_t : instance of Transform |
| 322 | New ``dev_head_t`` transformation using the averaged good head positions. |
| 323 | |
| 324 | Notes |
| 325 | ----- |
| 326 | .. versionchanged:: 1.7 |
| 327 | Support for multiple raw instances and position arrays was added. |
| 328 | """ |
| 329 | # Get weighted head pos trans and rot |
| 330 | if not isinstance(raw, list | tuple): |
| 331 | raw = [raw] |
| 332 | if not isinstance(pos, list | tuple): |
| 333 | pos = [pos] |
| 334 | if len(pos) != len(raw): |
| 335 | raise ValueError( |
| 336 | f"Number of head positions ({len(pos)}) must match the number of raw " |
| 337 | f"instances ({len(raw)})" |
| 338 | ) |
| 339 | hp = list() |
| 340 | dt = list() |
| 341 | for ri, (r, p) in enumerate(zip(raw, pos)): |
| 342 | _validate_type(r, BaseRaw, f"raw[{ri}]") |
| 343 | _validate_type(p, np.ndarray, f"pos[{ri}]") |
| 344 | hp_, dt_ = _raw_hp_weights(r, p) |
| 345 | hp.append(hp_) |
| 346 | dt.append(dt_) |
| 347 | hp = np.concatenate(hp, axis=0) |
| 348 | dt = np.concatenate(dt, axis=0) |
| 349 | dt /= dt.sum() |
| 350 | best_q = _average_quats(hp[:, 1:4], weights=dt) |
| 351 | trans = np.eye(4) |
| 352 | trans[:3, :3] = quat_to_rot(best_q) |
| 353 | trans[:3, 3] = dt @ hp[:, 4:7] |
| 354 | dist = np.linalg.norm(trans[:3, 3]) |
| 355 | if dist > 1: # less than 1 meter is sane |
| 356 | warn(f"Implausible head position detected: {dist} meters from device origin") |
| 357 | dev_head_t = Transform("meg", "head", trans) |
| 358 | return dev_head_t |
| 359 | |
| 360 |