[BUG] Fix mne.viz.Brain.add_volume_labels matrix ordering bug (#11730)

Co-authored-by: Eric Larson <[email protected]>

doc/changes/latest.inc

@@ -29,6 +29,7 @@ Enhancements
 Bugs
 ~~~~
 - Extended test to highlight bug in :func:`mne.stats.permutation_t_test` (:gh:`11575` by :newcontrib:`Joshua Calder-Travis`)
+- Fix bug where :meth:`mne.viz.Brain.add_volume_labels` used an incorrect orientation (:gh:`11730` by `Alex Rockhill`_)
 - Fix bug with :func:`mne.forward.restrict_forward_to_label` where cortical patch information was not adjusted (:gh:`11694` by `Eric Larson`_)
 - Fix bug with PySide6 compatibility (:gh:`11721` by `Eric Larson`_)
 - Fix hanging interpreter with matplotlib figures using ``mne/viz/_mpl_figure.py`` in spyder console and jupyter notebooks`(:gh:`11696` by `Mathieu Scheltienne`_)

mne/surface.py

@@ -1917,20 +1917,23 @@ def _marching_cubes(image, level, smooth=0, fill_hole_size=None, use_flying_edge
             f"{level.size} elements"
         )
 
+    # vtkImageData indexes as slice, row, col (Z, Y, X):
+    # https://discourse.vtk.org/t/very-confused-about-imdata-matrix-index-order/6608/2
+    # We can accomplish this by raveling with order='F' later, so we might as
+    # well make a copy with Fortran order now.
+    # We also use double as passing integer types directly can be problematic!
+    image = np.array(image, dtype=float, order="F")
+    image_shape = image.shape
+
     # fill holes
     if fill_hole_size is not None:
-        image = image.copy()  # don't modify original
         for val in level:
             bin_image = image == val
             mask = image == 0  # don't go into other areas
             bin_image = binary_dilation(bin_image, iterations=fill_hole_size, mask=mask)
             image[bin_image] = val
 
-    # force double as passing integer types directly can be problematic!
-    image_shape = image.shape
-    # use order='A' to automatically detect when Fortran ordering is needed
-    data_vtk = numpy_to_vtk(image.ravel(order="A").astype(float), deep=True)
-    del image
+    data_vtk = numpy_to_vtk(image.ravel(order="F"), deep=False)
 
     mc = vtkDiscreteFlyingEdges3D() if use_flying_edges else vtkDiscreteMarchingCubes()
     # create image
@@ -1978,8 +1981,8 @@ def _marching_cubes(image, level, smooth=0, fill_hole_size=None, use_flying_edge
         polydata = geometry.GetOutput()
         rr = vtk_to_numpy(polydata.GetPoints().GetData())
         tris = vtk_to_numpy(polydata.GetPolys().GetConnectivityArray()).reshape(-1, 3)
-        rr = np.ascontiguousarray(rr[:, ::-1])
-        tris = np.ascontiguousarray(tris[:, ::-1])
+        rr = np.ascontiguousarray(rr)
+        tris = np.ascontiguousarray(tris)
         out.append((rr, tris))
     return out
 

mne/tests/test_surface.py

@@ -245,12 +245,13 @@ def test_normal_orth():
 
 # 0.06 s locally even with all these params
 @pytest.mark.parametrize("dtype", (np.float64, np.uint16, ">i4"))
+@pytest.mark.parametrize("order", "FC")
 @pytest.mark.parametrize("value", (1, 12))
 @pytest.mark.parametrize("smooth", (0, 0.9))
-def test_marching_cubes(dtype, value, smooth):
+def test_marching_cubes(dtype, value, smooth, order):
     """Test creating surfaces via marching cubes."""
     pytest.importorskip("pyvista")
-    data = np.zeros((50, 50, 50), dtype=dtype)
+    data = np.zeros((50, 50, 50), dtype=dtype, order=order)
     data[20:30, 20:30, 20:30] = value
     level = [value]
     out = _marching_cubes(data, level, smooth=smooth)
@@ -260,8 +261,7 @@ def test_marching_cubes(dtype, value, smooth):
     rtol = 1e-2 if smooth else 1e-9
     assert_allclose(verts.sum(axis=0), [14700, 14700, 14700], rtol=rtol)
     tri_sum = triangles.sum(axis=0).tolist()
-    # old VTK (9.2.6), new VTK
-    assert tri_sum in [[363402, 360865, 350588], [364089, 359867, 350408]]
+    assert tri_sum in ([350588, 360865, 363402], [350408, 359867, 364089])
     # test fill holes
     data[24:27, 24:27, 24:27] = 0
     verts, triangles = _marching_cubes(data, level, smooth=smooth, fill_hole_size=2)[0]