`compute_average_dev_head_t` averages head rotation with squared time weights (translation uses linear weights)

Abierto Apto para principiantes
  • ¿Algo no está bien?
#14,478 2 comentarios 0 reacciones 0 asignados Ver en GitHub

Los mantenedores suelen responder en 1 día

Nadie ha tomado este issue todavía.

Evaluación

Dificultad
2/5
Tiempo estimado
1-3 horas
Aptitud para principiantes
66/100
Tipo de issue
Error
Claridad
Bien especificado
Estado de actividad
Activo
Stack tecnológico
numpy, python
Área
data

Línea de trabajo

El error está en _average_quats en mne/transforms.py, que es llamada por compute_average_dev_head_t en mne/preprocessing. Empieza por las líneas que multiplican los cuaterniones por los pesos antes del einsum que construye A, y comprueba dónde se elevan los pesos al cuadrado. Ejecuta el script de reproducción del issue: el ángulo de rotación debería salir en unos 0,1501 rad, no 0,1801. Está terminado cuando el ángulo de rotación coincide con la media ponderada y las pruebas existentes de transformaciones siguen pasando.

Escrito por el modelo de indexación a partir del texto del issue.

Descripción

BUG
Description of the problem

_average_quats (used by mne.preprocessing.compute_average_dev_head_t) is documented in a comment as the weighted quaternion average of Markley et al., which maximizes q^T A q with A = sum_i w_i q_i q_i^T. The implementation multiplies each quaternion by its weight and then forms the outer products, so it actually builds A = sum_i w_i^2 q_i q_i^T. The rotation is therefore averaged with squared weights, while the translation in the same function is averaged with the plain time weights, so the two parts of the returned transform are inconsistent.

Steps to reproduce
import numpy as np, mne
from mne.transforms import rot_to_quat

def rz(t):
    return np.array([[np.cos(t), -np.sin(t), 0], [np.sin(t), np.cos(t), 0], [0, 0, 1]])

sfreq = 100.0
raw = mne.io.RawArray(np.zeros((1, 1000)), mne.create_info(1, sfreq, "eeg"))
q = rot_to_quat(np.array([rz(0.0), rz(0.2)]))
# columns: time, q1, q2, q3, x, y, z, g, error, velocity
pos = np.array([[0.0, *q[0], 0.00, 0, 0, 1, 0, 0],
                [2.5, *q[1], 0.01, 0, 0, 1, 0, 0]])
trans = mne.preprocessing.compute_average_dev_head_t(raw, pos)["trans"]
print(np.arctan2(trans[1, 0], trans[0, 0]))  # 0.1801
print(trans[0, 3])                           # 0.0075
Link to data

No response

Expected results

The translation is the time-weighted mean (0.75 * 0.01 = 0.0075). The rotation should be the matching weighted quaternion mean: with A = sum_i w_i q_i q_i^T the dominant eigenvector corresponds to a rotation angle of 0.1501 rad (close to the weighted mean angle 0.75 * 0.2 = 0.15).

Actual results

The rotation angle is 0.1801 rad, which matches what squared weights give (0.5625 / (0.0625 + 0.5625) * 0.2 = 0.18). The head-position average is biased towards the longest segment.

Additional information

I suspect the problem is in mne/transforms.py, _average_quats:

quats = np.concatenate((_quat_real(quats)[..., np.newaxis], quats), -1)
quats *= weights[:, np.newaxis]
A = np.einsum("ij,ik->jk", quats, quats)  # sum of outer product of each q
Lenguaje dominante
Python
Estrellas
3.6k
Forks
1.6k
Merge medio
1 d 20 h
PR fusionados (30 d)
73

Preparar el entorno

Primeros pasos

  1. Lee el issue completo y luego la guía de contribución del proyecto.
  2. Comenta en el issue que vas a ocuparte — evita que dos personas hagan lo mismo.
  3. Haz un fork del repositorio y trabaja en una rama.
  4. Abre un pull request que haga referencia al número del issue.

Más de mne-tools/mne-python

Todos los issues de mne-tools/mne-python

Issues similares

Más issues de Python

Recibe los nuevos issues en tu correo

Un resumen breve de issues de GitHub para principiantes.