BUG: Fix bug with least-squares sphere fit (#13178)

Co-authored-by: autofix-ci[bot] <114827586+autofix-ci[bot]@users.noreply.github.com>

Author: larsoner

doc/changes/devel/13178.bugfix.rst

@@ -0,0 +1 @@
+Fix bug with least-squares fitting of head origin using digitization points in :func:`mne.preprocessing.maxwell_filter`, by `Eric Larson`_.

mne/_fiff/tests/test_reference.py

@@ -312,7 +312,11 @@ def test_set_eeg_reference_rest():
     same = [raw.ch_names.index(raw.info["bads"][0])]
     picks = np.setdiff1d(np.arange(len(raw.ch_names)), same)
     trans = None
-    sphere = make_sphere_model("auto", "auto", raw.info)
+    # Use fixed values from old sphere fit to reduce lines changed with fixed algorithm
+    sphere = make_sphere_model(
+        [-0.00413508, 0.01598787, 0.05175598],
+        0.09100286249131773,
+    )
     src = setup_volume_source_space(pos=20.0, sphere=sphere, exclude=30.0)
     assert src[0]["nuse"] == 223  # low but fast
     fwd = make_forward_solution(raw.info, trans, src, sphere)

mne/bem.py

@@ -1073,10 +1073,10 @@ def get_fitting_dig(info, dig_kinds="auto", exclude_frontal=True, verbose=None):
 
 
 @verbose
-def _fit_sphere_to_headshape(info, dig_kinds, verbose=None):
+def _fit_sphere_to_headshape(info, dig_kinds, *, verbose=None):
     """Fit a sphere to the given head shape."""
     hsp = get_fitting_dig(info, dig_kinds)
-    radius, origin_head = _fit_sphere(np.array(hsp), disp=False)
+    radius, origin_head = _fit_sphere(np.array(hsp))
     # compute origin in device coordinates
     dev_head_t = info["dev_head_t"]
     if dev_head_t is None:
@@ -1105,36 +1105,16 @@ def _fit_sphere_to_headshape(info, dig_kinds, verbose=None):
     return radius, origin_head, origin_device
 
 
-def _fit_sphere(points, disp="auto"):
+def _fit_sphere(points):
     """Fit a sphere to an arbitrary set of points."""
-    if isinstance(disp, str) and disp == "auto":
-        disp = True if logger.level <= 20 else False
-    # initial guess for center and radius
-    radii = (np.max(points, axis=1) - np.min(points, axis=1)) / 2.0
-    radius_init = radii.mean()
-    center_init = np.median(points, axis=0)
-
-    # optimization
-    x0 = np.concatenate([center_init, [radius_init]])
-
-    def cost_fun(center_rad):
-        d = np.linalg.norm(points - center_rad[:3], axis=1) - center_rad[3]
-        d *= d
-        return d.sum()
-
-    def constraint(center_rad):
-        return center_rad[3]  # radius must be >= 0
-
-    x_opt = fmin_cobyla(
-        cost_fun,
-        x0,
-        constraint,
-        rhobeg=radius_init,
-        rhoend=radius_init * 1e-6,
-        disp=disp,
-    )
-
-    origin, radius = x_opt[:3], x_opt[3]
+    # linear least-squares sphere fit, see for example
+    # https://stackoverflow.com/a/78909044
+    # TODO: At some point we should maybe reject outliers first...
+    A = np.c_[2 * points, np.ones((len(points), 1))]
+    b = (points**2).sum(axis=1)
+    x, _, _, _ = np.linalg.lstsq(A, b, rcond=1e-6)
+    origin = x[:3]
+    radius = np.sqrt(x[0] ** 2 + x[1] ** 2 + x[2] ** 2 + x[3])
     return radius, origin
 
 

mne/channels/tests/test_montage.py

@@ -1453,7 +1453,7 @@ def renamer(x):
     assert (orig_pos != new_pos).all()
 
     r0 = _fit_sphere(new_pos)[1]
-    assert_allclose(r0, [-0.001021, 0.014554, 0.041404], atol=1e-4)
+    assert_allclose(r0, [-0.001043, 0.01469, 0.041448], atol=1e-4)
     # spot check: Fp1 and Fpz
     assert_allclose(
         new_pos[:2],

mne/dipole.py

@@ -1573,7 +1573,7 @@ def fit_dipole(
         # Find the best-fitting sphere
         inner_skull = _bem_find_surface(bem, "inner_skull")
         inner_skull = inner_skull.copy()
-        R, r0 = _fit_sphere(inner_skull["rr"], disp=False)
+        R, r0 = _fit_sphere(inner_skull["rr"])
         # r0 back to head frame for logging
         r0 = apply_trans(mri_head_t["trans"], r0[np.newaxis, :])[0]
         inner_skull["r0"] = r0

mne/io/curry/tests/test_curry.py

@@ -280,7 +280,7 @@ def test_read_raw_curry_rfDC(fname, tol, mock_dev_head_t, tmp_path):
     assert n_other == 15
     pos = np.array(pos)
     nn = np.array(nn)
-    rad, origin = _fit_sphere(pos, disp=False)
+    rad, origin = _fit_sphere(pos)
     assert 0.11 < rad < 0.13
     pos -= origin
     pos /= np.linalg.norm(pos, axis=1, keepdims=True)

mne/preprocessing/tests/test_fine_cal.py

@@ -243,7 +243,7 @@ def test_fine_cal_systems(system, tmp_path):
         int_order = 5
         corrs = (0.13, 0.0, 0.12)
         sfs = [4, 5, 125, 155]
-        corr_tol = 0.3
+        corr_tol = 0.34
     else:
         assert system == "triux", f"Unknown system {system}"
         raw = read_raw_fif(tri_fname)

mne/preprocessing/tests/test_maxwell.py

@@ -173,6 +173,7 @@
 
 def _assert_n_free(raw_sss, lower, upper=None):
     """Check the DOF."""
+    __tracebackhide__ = True
     upper = lower if upper is None else upper
     n_free = raw_sss.info["proc_history"][0]["max_info"]["sss_info"]["nfree"]
     assert lower <= n_free <= upper, f"nfree fail: {lower} <= {n_free} <= {upper}"
@@ -367,10 +368,9 @@ def test_other_systems():
     hsp_path = kit_dir / "test_hsp.txt"
     raw_kit = read_raw_kit(sqd_path, str(mrk_path), str(elp_path), str(hsp_path))
     with (
-        pytest.warns(RuntimeWarning, match="more than 20 mm from head frame origin"),
         pytest.raises(NotImplementedError, match="Cannot create forward solution with"),
     ):
-        maxwell_filter(raw_kit)
+        maxwell_filter(raw_kit, verbose=True)
     with catch_logging() as log:
         raw_sss = maxwell_filter(
             raw_kit, origin=(0.0, 0.0, 0.04), ignore_ref=True, verbose=True
@@ -381,24 +381,21 @@ def test_other_systems():
         raw_kit, origin=(0.0, 0.0, 0.04), ignore_ref=True, mag_scale="auto"
     )
     assert_allclose(raw_sss._data, raw_sss_auto._data)
-    # The KIT origin fit is terrible
-    with pytest.warns(RuntimeWarning, match="more than 20 mm"):
-        with catch_logging() as log:
-            pytest.raises(
-                RuntimeError, maxwell_filter, raw_kit, ignore_ref=True, regularize=None
-            )  # bad condition
-            raw_sss = maxwell_filter(
-                raw_kit,
-                origin="auto",
-                ignore_ref=True,
-                bad_condition="info",
-                verbose=True,
-            )
+    with catch_logging() as log:
+        pytest.raises(
+            RuntimeError, maxwell_filter, raw_kit, ignore_ref=True, regularize=None
+        )  # bad condition
+        raw_sss = maxwell_filter(
+            raw_kit,
+            origin="auto",
+            ignore_ref=True,
+            bad_condition="info",
+            verbose=True,
+        )
     log = log.getvalue()
-    assert "badly conditioned" in log
-    assert "more than 20 mm from" in log
-    # fits can differ slightly based on scipy version, so be lenient here
-    _assert_n_free(raw_sss, 28, 34)  # bad origin == brutal reg
+    assert "badly conditioned" not in log
+    assert "more than 20 mm from" not in log
+    _assert_n_free(raw_sss, 67, 67)
     # Let's set the origin
     with catch_logging() as log:
         raw_sss = maxwell_filter(

mne/simulation/tests/test_raw.py

@@ -220,7 +220,11 @@ def raw_data():
 
     src = read_source_spaces(src_fname)
     trans = read_trans(trans_fname)
-    sphere = make_sphere_model("auto", "auto", raw.info)
+    # Use fixed values from old sphere fit to reduce lines changed with fixed algorithm
+    sphere = make_sphere_model(
+        [-0.00413508, 0.01598787, 0.05175598],
+        0.09100286249131773,
+    )
     stc = _make_stc(raw, src)
     return raw, src, stc, trans, sphere
 
@@ -385,9 +389,9 @@ def test_simulate_raw_bem(raw_data):
     locs = np.concatenate([s["rr"][s["vertno"]] for s in src])
     tmax = (len(locs) - 1) / raw.info["sfreq"]
     cov = make_ad_hoc_cov(raw.info)
-    # The tolerance for the BEM is surprisingly high (28) but I get the same
+    # The tolerance for the BEM is surprisingly high but I get the same
     # result when using MNE-C and Xfit, even when using a proper 5120 BEM :(
-    for use_raw, bem, tol in ((raw_sim_sph, sphere, 3), (raw_sim_bem, bem_fname, 31)):
+    for use_raw, bem, tol in ((raw_sim_sph, sphere, 4), (raw_sim_bem, bem_fname, 31)):
         events = find_events(use_raw, "STI 014")
         assert len(locs) == 6
         evoked = Epochs(use_raw, events, 1, 0, tmax, baseline=None).average()

mne/surface.py

@@ -220,7 +220,7 @@ def get_meg_helmet_surf(info, trans=None, *, verbose=None):
         )
         hull = ConvexHull(rr)
         rr = rr[np.unique(hull.simplices)]
-        R, center = _fit_sphere(rr, disp=False)
+        R, center = _fit_sphere(rr)
         sph = _cart_to_sph(rr - center)[:, 1:]
         # add a point at the front of the helmet (where the face should be):
         # 90 deg az and maximal el (down from Z/up axis)