Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 10 commits into from
Jun 14, 2023

Conversation

alexrockhill
Copy link
Contributor

Fixes #11728

@alexrockhill
Copy link
Contributor Author

This works for me now

import os.path as op
import mne
mne.viz.set_3d_backend('pyvistaqt') # GUI requires pyvista backend
mne.viz.set_3d_options(depth_peeling=False, antialias=False, multi_samples=1) # for bugs of pure black plot of 3D viz

fsaverage_dir = op.dirname(mne.datasets.fetch_fsaverage())
mne.datasets.fetch_hcp_mmp_parcellation()

brain = mne.viz.Brain(
    "fsaverage",
    subjects_dir=fsaverage_dir,
    cortex='low_contrast', alpha=0.2, background='white'
)

brain.add_volume_labels(aseg="aparc+aseg")

Comment on lines 2647 to 2654
ornt = nib.orientations.axcodes2ornt(
nib.orientations.aff2axcodes(aseg.affine)
).astype(int)
ras_ornt = nib.orientations.axcodes2ornt("RAS")
ornt_trans = nib.orientations.ornt_transform(ornt, ras_ornt)
aseg_data = nib.orientations.apply_orientation(aseg_data, ornt_trans)
aff_trans = nib.orientations.inv_ornt_aff(ornt_trans, aseg.shape)
vox_mri_t = np.dot(vox_mri_t, aff_trans)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am confused why this is necessary -- the vox_mri_t should already take care of putting everything in MRI surface RAS coordinates for us. Also, there are a number of places we use vox_mri_t and don't use any axcodes2ornt stuff, so it's weird this is the one place we need to do it.

Copy link
Member

@larsoner larsoner Jun 12, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

... to me it suggests maybe there is some bug with how _marching_cubes works or something, rather than some adjustment we need to make to vox_mri_t.

I would investigate deeper by using freeview and some specific voxel i, j, k in some known subvolume and what their viewer says for MRI surface RAS x, y, z, confirm that this is just vox_mri_t @ (i, j, k, 1) (it had better be!), then inspect the outputs at various points of this function for that label. It should tell us where the problem is...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So it looks like it was caused by the fortran order check in _marching_cubes. I think the solution I just pushed works, here is one issue where the order needed to be fortran #10860 which lead to the order being auto-detected with "A". Checking to be sure that that works now too.

@alexrockhill
Copy link
Contributor Author

Huh so oddly, in the tumor case, the np.array.flags.c_contiguous is True and f_contiguous is False but order="F" makes the tumor look normal and in the case of this issue everything is opposite. That would suggest something like order="F" if image.flags.c_contiguous else "C" which looks very confusing but I think would work.

@alexrockhill
Copy link
Contributor Author

alexrockhill commented Jun 12, 2023

import os.path as op
import mne
mne.viz.set_3d_backend('pyvistaqt') # GUI requires pyvista backend
mne.viz.set_3d_options(depth_peeling=False, antialias=False, multi_samples=1) # for bugs of pure black plot of 3D viz
 
fsaverage_dir = op.dirname(mne.datasets.fetch_fsaverage())
mne.datasets.fetch_hcp_mmp_parcellation()
  
brain = mne.viz.Brain(
     "fsaverage",
      subjects_dir=fsaverage_dir,
      cortex='low_contrast', alpha=0.2, background='white'
) 
brain.add_volume_labels(aseg="aparc+aseg")
import nrrd
import mne
from mne.surface import _marching_cubes

data, header = nrrd.read('tumor.nrrd', index_order='F')
verts, faces = _marching_cubes(data, [1], smooth=0)[0]
fig = mne.viz.create_3d_figure((600, 600), scene=False)
fig.mesh(*verts.T, faces, 'red')
fig.show()

Both of these work now using the tumor.nrrd file from here: #10860. It's not ideal that it's opposite of what you'd expect but I'm not really sure if there's a workaround that's better other than just documenting it with a verbose comment.

@alexrockhill
Copy link
Contributor Author

Okay, after much debugging, part of the issue was with image = image.copy() which made a shallow copy that changed the index ordering from Fortran to C without actually changing the index ordering. I checked and this version works for the above two MWE as well as the sEEG example and https://mne.tools/dev/auto_tutorials/inverse/60_visualize_stc.html#sphx-glr-auto-tutorials-inverse-60-visualize-stc-py which are the two uses of add_volume_labels and the latter looks like it's currently broken on main so this will fix that. Should be good to go.

@alexrockhill alexrockhill changed the title [BUG] Fix add volume labels, reorient to ras [BUG] Fix mne.viz.Brain.add_volume_labels matrix ordering bug Jun 12, 2023
@alexrockhill
Copy link
Contributor Author

Oh also I in doing this I read into the vtk functions and it says that if you delete the reference image in a numpy_to_vtk you can cause a segfault so I didn't do that. We might do that other places so maybe a good FYI @larsoner.

numpy_to_vtk(num_array, deep=0, array_type=None)
    Converts a real numpy Array to a VTK array object.
    
    This function only works for real arrays.
    Complex arrays are NOT handled.  It also works for multi-component
    arrays.  However, only 1, and 2 dimensional arrays are supported.
    This function is very efficient, so large arrays should not be a
    problem.
    
    If the second argument is set to 1, the array is deep-copied from
    from numpy. This is not as efficient as the default behavior
    (shallow copy) and uses more memory but detaches the two arrays
    such that the numpy array can be released.
    
    WARNING: You must maintain a reference to the passed numpy array, if
    the numpy data is gc'd and VTK will point to garbage which will in
    the best case give you a segfault.
    
    Parameters:
    
    num_array
      a 1D or 2D, real numpy array.

@larsoner larsoner enabled auto-merge (squash) June 13, 2023 23:52
@larsoner larsoner merged commit 6fd91f6 into mne-tools:main Jun 14, 2023
@larsoner
Copy link
Member

Thanks @alexrockhill !

@Dr-Chen-Xiaoyu
Copy link

The lastest dev version of MNE correctly plots. 👍👍👍 Thanks a lot !!!

@alexrockhill alexrockhill deleted the cf branch June 15, 2023 20:03
larsoner added a commit to larsoner/mne-python that referenced this pull request Jun 23, 2023
* upstream/main: (24 commits)
  Allow int-like as ID of make_fixed_length_events (mne-tools#11748)
  Easycap-M43 montage (mne-tools#11744)
  ENH: Create a Calibrations class for eyetracking data (mne-tools#11719)
  Fix alphabetical order in overview/people.rst, fix sphinx formatting in docstrings and set verbose to keyword-only (mne-tools#11745)
  Add Mathieu Scheltienne to MNE-Python Steering Council (mne-tools#11741)
  removed requirement for curv.*h files to create Brain object (mne-tools#11704)
  [BUG] Fix mne.viz.Brain.add_volume_labels matrix ordering bug (mne-tools#11730)
  Fix installer links (mne-tools#11729)
  MAINT: Update for PyVista deprecation (mne-tools#11727)
  MAINT: Update roadmap (mne-tools#11724)
  MAINT: Update download link [skip azp] [skip cirrus] [skip actions]
  fix case for chpi_info[1] == None (mne-tools#11714)
  Add cmap argument for mne.viz.utils.plot_sensors (mne-tools#11720)
  BUG: Fix one more PySide6 bug (mne-tools#11723)
  MAINT: Fix PySide6 and PyVista compat (mne-tools#11721)
  MRG: If _check_fname() cannot find a file, display the path in quotation marks to help spot accidental trailing spaces (mne-tools#11718)
  Add "array-like" to `_validate_type()` (mne-tools#11713)
  MAINT: Avoid problematic PySide6 (mne-tools#11715)
  Fix installer links (mne-tools#11709)
  Updating change log after PR mne-tools#11575 (mne-tools#11707)
  ...
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Issue in brain.add_volume_labels()
3 participants