nipy / PySurfer

Cortical neuroimaging visualization in Python
https://pysurfer.github.io/
BSD 3-Clause "New" or "Revised" License
240 stars 97 forks source link

[MRG]: Set surface #191

Closed christianbrodbeck closed 7 years ago

christianbrodbeck commented 7 years ago

I often find myself plotting with one surface and then wanting to look at it with another. This implementation is based on keeping track of all the meshes that are added to the brain. I've implemented it for data, labels and annotations, not yet for all other layer types because I wanted to confirm first that there is no better way to implement it. I was wondering whether all the mesh objects could be combined into one somehow (re #177), but it looks like they all need to have a different scalars attribute, and I don't see a straightforward way to separate the scalars from the triangular mesh...

mwaskom commented 7 years ago

Cool!

My first thought is that it should probably be set_surf to be more consistent with the keyword in the Brain constructor (and the terminology in most of Freesurfer.

coveralls commented 7 years ago

Coverage Status

Coverage increased (+0.09%) to 78.509% when pulling 9321b587570ba787a0cece6fcd7d4ab29a62b82d on christianbrodbeck:set_surface into f94bea215be1605b6bbea2b0b98af772d30c0a62 on nipy:master.

coveralls commented 7 years ago

Coverage Status

Coverage increased (+0.09%) to 78.509% when pulling 188732ec584a1b741d2f38b05b39e6ea855351ab on christianbrodbeck:set_surface into f94bea215be1605b6bbea2b0b98af772d30c0a62 on nipy:master.

christianbrodbeck commented 7 years ago

I'm working on sharing data based on this. I ran into an error with Threshold filters which try to access a nonexistent .point_data attribute:

screen shot 2017-06-13 at 8 12 34 pm

Then I noticed a similar error happens when I I run the second (longer) example from the link. Is this a bug in Mayavi 4.5? Can others run the example without getting error dialogs?

I also had to set the backend for many tests to auto to avoid getting this dialog:

screen shot 2017-06-13 at 8 25 39 pm

and this error:

apptools.preferences.preferences: DEBUG: loading preferences from </Users/christian/.enthought/mayavi_e3/preferences.ini>
apptools.preferences.preferences: DEBUG: loading preferences from <<open file '//anaconda/envs/mayavi45/lib/python2.7/site-packages/mayavi/preferences/preferences.ini', mode 'rb' at 0x1188470c0>>
apptools.preferences.preferences: DEBUG: loading preferences from <<open file '//anaconda/envs/mayavi45/lib/python2.7/site-packages/tvtk/plugins/scene/preferences.ini', mode 'rb' at 0x1188478a0>>
mayavi.core.registry: DEBUG: Engine [<mayavi.core.null_engine.NullEngine object at 0x102d38e30>] named NullEngine1 registered
mayavi.core.common: ERROR: No valid current object, please select an active object.
traits: ERROR: Exception occurred in traits notification handler for object: <mayavi.tools.filters.SetActiveAttributeFactory object at 0x1029d8b30>, trait: point_scalars, old value: , new value: 2
Traceback (most recent call last):
  File "//anaconda/envs/mayavi45/lib/python2.7/site-packages/traits/trait_notifiers.py", line 340, in __call__
    self.handler( *args )
  File "//anaconda/envs/mayavi45/lib/python2.7/site-packages/mayavi/tools/pipe_base.py", line 198, in _anytrait_changed
    setattr(obj, components[-1], value)
  File "//anaconda/envs/mayavi45/lib/python2.7/site-packages/mayavi/core/trait_defs.py", line 97, in set_value
    raise TraitError(object, name,"one of %s"%values, value)
TraitError: The 'point_scalars_name' trait of a SetActiveAttribute instance must be one of [], but a value of '2' <type 'str'> was specified.
christianbrodbeck commented 7 years ago

I've opened an issue at https://github.com/enthought/mayavi/issues/511

larsoner commented 7 years ago

I think 4.4.4 (or whatever it was) should still live in the Anaconda ecosystem somewhere, can you see if you can run it on older Mayavi version(s)?

That blog post is over 7 years old, so I'm not too surprised it no longer works :(

larsoner commented 7 years ago

On my OSX machine I get no errors. On my Linux machine, I get: screenshot from 2017-06-14 10-44-40

Followed by: screenshot from 2017-06-14 10-44-56

with messages:

No handlers could be found for logger "mayavi.core.common"
Exception occurred in traits notification handler.
Please check the log file for details.
Exception occurred in traits notification handler for object: <mayavi.tools.decorations.Axes object at 0x7f6c2b295c50>, trait: extent, old value: [0 0 0 0 0 0], new value: [ -50.   49.  -50.   49.    0.  100.]
Traceback (most recent call last):
  File "/home/larsoner/anaconda2/lib/python2.7/site-packages/traits/trait_notifiers.py", line 340, in __call__
    self.handler( *args )
  File "/home/larsoner/anaconda2/lib/python2.7/site-packages/mayavi/tools/decorations.py", line 373, in _extent_changed
    axes.module_manager.source.outputs[0].bounds
AttributeError: 'AssignAttribute' object has no attribute 'bounds'

But the plot works (here I've replaced w = np.random.RandomState(0).randn(*z.shape)): screenshot from 2017-06-14 10-44-07

And I get no errors if I remove the mlab.axes() call. Can you see if it's the case for you?

christianbrodbeck commented 7 years ago

Yes you're right, it was just the mlab.axes() call!

So now for the actual problem: If I replace the surf = ... line with

thresh = mlab.pipeline.threshold(active_attr, low=0.5)
surf = mlab.pipeline.surface(thresh)

I get:

screen shot 2017-06-14 at 11 07 25 am

then

screen shot 2017-06-14 at 11 07 28 am

and

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
/Users/christian/Documents/scripts/mayavi/scalar_data.py in <module>()
     34 ###
     35 thresh = mlab.pipeline.threshold(active_attr, low=0.5)
---> 36 surf = mlab.pipeline.surface(thresh)
     37 ###
     38 

/Users/christian/anaconda/envs/mayavi/lib/python2.7/site-packages/mayavi/tools/pipe_base.pyc in the_function(*args, **kwargs)
     36 def make_function(factory_class):
     37     def the_function(*args, **kwargs):
---> 38         factory = factory_class(*args, **kwargs)
     39         return factory._target
     40 

/Users/christian/anaconda/envs/mayavi/lib/python2.7/site-packages/mayavi/tools/modules.pyc in __init__(self, *args, **kwargs)
    154 
    155     def __init__(self, *args, **kwargs):
--> 156         super(DataModuleFactory, self).__init__(*args, **kwargs)
    157         # We are adding data to the scene, reset the zoom:
    158         scene = self._scene.scene

/Users/christian/anaconda/envs/mayavi/lib/python2.7/site-packages/mayavi/tools/pipe_base.pyc in __init__(self, parent, **kwargs)
    161         # Now calling the traits setter, so that traits handlers are
    162         # called
--> 163         self.set(**traits)
    164         if scene is not None:
    165             scene.disable_render = not self._do_redraw

/Users/christian/anaconda/envs/mayavi/lib/python2.7/site-packages/mayavi/tools/pipe_base.pyc in set(self, trait_change_notify, **traits)
    177             try:
    178                 if callback is not None:
--> 179                     callback()
    180                 self._anytrait_changed(trait, value)
    181             except TraitError:

/Users/christian/anaconda/envs/mayavi/lib/python2.7/site-packages/mayavi/tools/modules.pyc in _colormap_changed(self)
    123             self._target.module_manager.scalar_lut_manager.reverse_lut = True
    124             self._target.module_manager.vector_lut_manager.reverse_lut = True
--> 125         self._target.module_manager.scalar_lut_manager.lut_mode = colormap
    126         self._target.module_manager.vector_lut_manager.lut_mode = colormap
    127 

AttributeError: 'NoneType' object has no attribute 'scalar_lut_manager'
larsoner commented 7 years ago

I see the same error. Can you confirm that this fixes it? Change mayavi/filters/threshold.py line ~223:

    def _get_data_range(self):
        """Returns the range of the input scalar data."""
        input = self.inputs[0].outputs[0]
        data_range = []
        if hasattr(input, 'point_data'):
            ps = input.point_data.scalars
        else:
            ps = None
        if hasattr(input, 'cell_data'):
            cs = input.cell_data.scalars
        else:
            cs = None
larsoner commented 7 years ago

(The Mayavi people have a ton of issues open, but you could open a PR as a potential fix, including some minimal example -- it makes it more likely for a proper fix to land)

christianbrodbeck commented 7 years ago

With the dev version of Mayavi I don't get this error, so I assume it has been fixed already (adding the missing attributes to the other filters)...

So for this PR, I assume the best course is to build in a workaround with a comment about what can be changed once mayavi updates?

larsoner commented 7 years ago

So for this PR, I assume the best course is to build in a workaround with a comment about what can be changed once mayavi updates?

Can we work around it by checking the mayavi version and monkey-patching a None?

christianbrodbeck commented 7 years ago

While updating all the displays I noticed that overlays are handled in a way that is not quite consistent with other objects, by being implemented with two classes (one to represent data, one to display it) while all other display types are implemented with dictionaries that are managed by Brain and _Hemisphere. This also seems to lead to inconsistent API in which annotations are removed through Brain.remove_annot() while overlays are removed with Brain.overlays_dict[x].remove(). Is specialized classes for displays something we want to move towards or way from? While it would be easier to reduce reliance on the class (keep it for backwards compatibility or deprecate it) I think classes could make things somewhat more transparent, and possibly expose functions that make it easier for users to customize displays.

larsoner commented 7 years ago

all other display types are implemented with dictionaries that are managed by Brain and _Hemisphere

This sounds like it's more backend-agnostic, which is something that I'd like to promote. Think about swapping in e.g. a matplotlib rendering backend, it would be easier to do if things are managed by Brain and _Hemisphere right?

annotations are removed through Brain.remove_annot() while overlays are removed with Brain.overlays_dict[x].remove()

Perhaps this is because people only ever thought about having a single annot at a time, whereas multiple overlays is a reasonable thing to use?

Is specialized classes for displays something we want to move towards or way from?

Can you write a few lines of (1) the current API and (2) the proposed alternatives APIs that move toward/away from the specialized classes? I still can't quite see it.

mwaskom commented 7 years ago

While updating all the displays I noticed that overlays are handled in a way that is not quite consistent with other objects, by being implemented with two classes (one to represent data, one to display it) while all other display types are implemented with dictionaries that are managed by Brain and _Hemisphere.

This is very, very old code. IIRC the motivation here was that overlays are a little weird in that you are mapping a single file to multiple mayavi piplines (as the positive and negative components are shown independently). That motivated creating a separate class for Overlays instead of just sticking the mayavi objects in a dictionary as is otherwise the case.

Perhaps this is because people only ever thought about having a single annot at a time, whereas multiple overlays is a reasonable thing to use?

This is correct though less directly relevant to having a class define the display type — i.e. you can have multiple outputs from add_data but they are just stored in a list I think.

christianbrodbeck commented 7 years ago

After updating the overlay object I think this PR should not be affected the classes vs. dicts discussion, so if you want to keep it open I would move it into a separate issue; but here my thoughts:

IIRC the motivation here was that overlays are a little weird in that you are mapping a single file to multiple mayavi piplines (as the positive and negative components are shown independently). That motivated creating a separate class for Overlays instead of just sticking the mayavi objects in a dictionary as is otherwise the case.

I guess these dictionaries grew over time in complexity... after adding a data-layer:

In [2]: brain.data_dict
Out[2]: 
{'lh': {'array': array([[ 1.43034017,  0.51767367,  0.45605102,  1.70507824],
         [ 0.26418859,  0.52732748,  0.90342617,  0.99659139],
         [ 1.78131807,  1.26297808,  1.32284141,  0.9324193 ],
         ..., 
         [ 2.65984082,  4.03173685,  6.1883111 ,  4.20004129],
         [ 2.9126246 ,  3.85630774,  3.45926285,  2.47084212],
         [ 3.17013526,  3.28790569,  2.09640741,  1.67996788]], dtype=float32),
  'colorbars': [<mayavi.core.lut_manager.LUTManager at 0x12b813950>],
  'fmax': 22.0,
  'fmid': 18.0,
  'fmin': 13.0,
  'layer': 0,
  'orig_ctable': array([[ 10,   0,   0, 255],
         [ 13,   0,   0, 255],
         [ 15,   0,   0, 255],
         ..., 
         [255, 255, 247, 255],
         [255, 255, 251, 255],
         [255, 255, 255, 255]], dtype=uint8),
  'smooth_mat': <163842x2562 sparse matrix of type '<type 'numpy.float64'>'
    with 847182 stored elements in COOrdinate format>,
  'smoothing_steps': 10,
  'surfaces': [<mayavi.modules.surface.Surface at 0x138464d10>],
  'time': array([ 0.08,  0.09,  0.1 ,  0.11]),
  'time_idx': 0,
  'time_label': <function __main__.<lambda>>,
  'transparent': True,
  'vertices': array([   0,    2,    3, ..., 1348, 2451, 2452], dtype=uint32)},
 'rh': {'array': array([[ 1.1138016 ,  1.63564575,  1.56338835,  0.62299478],
         [ 4.47758293,  4.73766136,  3.3510077 ,  1.50541949],
         [ 0.35411292,  0.7056011 ,  0.65926534,  0.89739203],
         ..., 
         [ 4.18813515,  3.81059504,  1.91492558,  1.60494661],
         [ 3.19541049,  4.48996067,  3.88375044,  2.12374163],
         [ 3.23499441,  6.16087055,  6.44084167,  3.88419032]], dtype=float32),
  'colorbars': [<mayavi.core.lut_manager.LUTManager at 0x12ef84fb0>],
  'fmax': 22.0,
  'fmid': 18.0,
  'fmin': 13.0,
  'layer': 1,
  'orig_ctable': array([[ 10,   0,   0, 255],
         [ 13,   0,   0, 255],
         [ 15,   0,   0, 255],
         ..., 
         [255, 255, 247, 255],
         [255, 255, 251, 255],
         [255, 255, 255, 255]], dtype=uint8),
  'smooth_mat': <163842x2562 sparse matrix of type '<type 'numpy.float64'>'
    with 847182 stored elements in COOrdinate format>,
  'smoothing_steps': 10,
  'surfaces': [<mayavi.modules.surface.Surface at 0x12ef84530>],
  'time': array([ 0.08,  0.09,  0.1 ,  0.11]),
  'time_idx': 0,
  'time_label': <function __main__.<lambda>>,
  'transparent': True,
  'vertices': array([   0,    2,    3, ..., 1348, 2451, 2452], dtype=uint32)}}

And this is only one layer... I assume dating back from when only one layer was possible. If there are multiple layers they are stored in Brain._data_dicts (which I made private when I added it to stop adding more attributes that have to be kept for backwards compatibility).

Can you write a few lines of (1) the current API and (2) the proposed alternatives APIs that move toward/away from the specialized classes? I still can't quite see it.

What I was thinking about (analogous to overlays) was replacing those dictionaries with a DataLayer object that could have methods like .scale_colormap(). As such I think it would actually move towards backend-neutrality because it would interpose a layer of PySurfer API in front of the mayavi objects' API rather than exposing the mayavi objects themselves.

larsoner commented 7 years ago

As such I think it would actually move towards backend-neutrality because it would interpose a layer of PySurfer API in front of the mayavi objects' API rather than exposing the mayavi objects themselves.

That sounds nice. But I agree it should be another PR if possible.

mwaskom commented 7 years ago

That sounds nice. But I agree it should be another PR if possible.

Agreed — it's a good idea but will be disruptive and will require a lot of thought to get right (which clearly wasn't done the first time, hence the mess).

coveralls commented 7 years ago

Coverage Status

Coverage increased (+0.1%) to 78.532% when pulling 3a09db16c8783111f998fbaa71d509c7b989e7b1 on christianbrodbeck:set_surface into f94bea215be1605b6bbea2b0b98af772d30c0a62 on nipy:master.

christianbrodbeck commented 7 years ago

Tests pass.

Regarding naming – I preferred .set_surface() because it's not an abbreviation, and in __init__ surf is a positional argument so the name is not highly salient... But if you prefer set_surf() I will update the name.

mwaskom commented 7 years ago

Regarding naming – I preferred .set_surface() because it's not an abbreviation, and in init surf is a positional argument so the name is not highly salient...

In my experience, the term surf is very commonly used in Freesurfer-land as a parameter that takes the an argument of white, pial, inflated, etc. that identifies different representations of the triangular mesh. Whereas "surface" is a more general term for the 2D space that is differentiated from the "volume". But others can weigh in.

christianbrodbeck commented 7 years ago

That makes sense, coming from the MNE side I did not know about this FreeSurfer convention...

christianbrodbeck commented 7 years ago

Search-and-replace also found some surface in the examples, I replaced those too for consistency.

coveralls commented 7 years ago

Coverage Status

Coverage increased (+0.1%) to 78.532% when pulling f2538348d88dbf3d1a54de1ab903a8da2fc45887 on christianbrodbeck:set_surface into f94bea215be1605b6bbea2b0b98af772d30c0a62 on nipy:master.

coveralls commented 7 years ago

Coverage Status

Coverage increased (+0.1%) to 78.532% when pulling 44bd47d45475f9099f4b26e327ea541211698f56 on christianbrodbeck:set_surface into f94bea215be1605b6bbea2b0b98af772d30c0a62 on nipy:master.

coveralls commented 7 years ago

Coverage Status

Coverage increased (+0.1%) to 78.532% when pulling e7d28a5b2db46b9f159f2e2e21e7788d5cc450d9 on christianbrodbeck:set_surface into f94bea215be1605b6bbea2b0b98af772d30c0a62 on nipy:master.

coveralls commented 7 years ago

Coverage Status

Coverage increased (+0.1%) to 78.532% when pulling bd99b944660310076b20d8741be8121ab612e2a2 on christianbrodbeck:set_surface into f94bea215be1605b6bbea2b0b98af772d30c0a62 on nipy:master.

coveralls commented 7 years ago

Coverage Status

Coverage increased (+0.1%) to 78.552% when pulling c88df561e9d81c125288ff7efb27a13ae86d3c3f on christianbrodbeck:set_surface into f94bea215be1605b6bbea2b0b98af772d30c0a62 on nipy:master.

christianbrodbeck commented 7 years ago

This is ready for review now. To test the new function use Brain.set_surf(...).

larsoner commented 7 years ago

I wonder if we should (via deprecation) explicitly set the focal point to 0, 0, 0 instead. It would change interaction a little bit but perhaps in a useful way.

christianbrodbeck commented 7 years ago

And align the geometry so that 0, 0, 0 is the actual center?

larsoner commented 7 years ago

Ideally we wouldn't need to do any realignment assuming the FreeSurfer surfaces are defined properly with 0, 0, 0 to have some meaning. But maybe they aren't, in which case aligning to the origin maybe isn't such a good idea.

christianbrodbeck commented 7 years ago

I think they aren't, there is a comment in the code under Brain.set_view() "blame"d on you :)

            # DO NOT set focal point, can screw up non-centered brains
            # view['focalpoint'] = (0.0, 0.0, 0.0)
larsoner commented 7 years ago

Whoops :)

Let's leave the resetting in for now, and maybe we can think of a better solution later

mwaskom commented 7 years ago

Good thinking, past Eric.

coveralls commented 7 years ago

Coverage Status

Coverage increased (+0.05%) to 78.479% when pulling da0986a8465cd3d8e2d70a431c3075c2329d6a38 on christianbrodbeck:set_surface into 9fa66d3f4313b4234a63adbd9ddaf51da3d76d50 on nipy:master.

coveralls commented 7 years ago

Coverage Status

Coverage increased (+0.05%) to 78.479% when pulling da0986a8465cd3d8e2d70a431c3075c2329d6a38 on christianbrodbeck:set_surface into 9fa66d3f4313b4234a63adbd9ddaf51da3d76d50 on nipy:master.

coveralls commented 7 years ago

Coverage Status

Coverage increased (+0.05%) to 78.479% when pulling da0986a8465cd3d8e2d70a431c3075c2329d6a38 on christianbrodbeck:set_surface into 9fa66d3f4313b4234a63adbd9ddaf51da3d76d50 on nipy:master.

coveralls commented 7 years ago

Coverage Status

Coverage increased (+0.05%) to 78.479% when pulling da0986a8465cd3d8e2d70a431c3075c2329d6a38 on christianbrodbeck:set_surface into 9fa66d3f4313b4234a63adbd9ddaf51da3d76d50 on nipy:master.

coveralls commented 7 years ago

Coverage Status

Coverage decreased (-24.4%) to 54.054% when pulling da0986a8465cd3d8e2d70a431c3075c2329d6a38 on christianbrodbeck:set_surface into 9fa66d3f4313b4234a63adbd9ddaf51da3d76d50 on nipy:master.

coveralls commented 7 years ago

Coverage Status

Coverage decreased (-24.4%) to 54.054% when pulling da0986a8465cd3d8e2d70a431c3075c2329d6a38 on christianbrodbeck:set_surface into 9fa66d3f4313b4234a63adbd9ddaf51da3d76d50 on nipy:master.

mwaskom commented 7 years ago

Uh oh; why did test coverage go down 24%?

christianbrodbeck commented 7 years ago

That seems wrong – I pushed too early and then force pushed a few times to amend, maybe that confused coverage – Travis shows 72% for this PR as well as for master.

christianbrodbeck commented 7 years ago

I successfully ran the tests with this env: $ conda create -n mayavi44 numpy scipy matplotlib mayavi=4.4 h5py pyqt=4.10 pygments=1 so it looks like monkey-patching for == 4.5 is right. I hope the new push will fix the coverage.

coveralls commented 7 years ago

Coverage Status

Coverage increased (+0.7%) to 79.088% when pulling 8c5222611685c2880c5d0c59fcd2e83643d38445 on christianbrodbeck:set_surface into 9fa66d3f4313b4234a63adbd9ddaf51da3d76d50 on nipy:master.

larsoner commented 7 years ago

This causes a freeze in MNE tests with "No valid current object" popup, can you try nosetests mne/viz/tests/test_3d.py and see if you can replicate?

larsoner commented 7 years ago

https://travis-ci.org/mne-tools/mne-python/jobs/246229590#L2360

larsoner commented 7 years ago
mayavi.core.common: ERROR: No valid current object, please select an active object.
traits: ERROR: Exception occurred in traits notification handler for object: <mayavi.tools.filters.SetActiveAttributeFactory object at 0x12f852ad0>, trait: point_scalars, old value: , new value: 2
Traceback (most recent call last):
  File "/Users/larsoner/anaconda/lib/python2.7/site-packages/traits/trait_notifiers.py", line 340, in __call__
    self.handler( *args )
  File "/Users/larsoner/anaconda/lib/python2.7/site-packages/mayavi/tools/pipe_base.py", line 198, in _anytrait_changed
    setattr(obj, components[-1], value)
  File "/Users/larsoner/anaconda/lib/python2.7/site-packages/mayavi/core/trait_defs.py", line 97, in set_value
    raise TraitError(object, name,"one of %s"%values, value)
TraitError: The 'point_scalars_name' trait of a SetActiveAttribute instance must be one of [], but a value of '2' <type 'str'> was specified.
christianbrodbeck commented 7 years ago

If I recall correctly that was the error I got with the test backend, the reason for changing them to auto. It seems to indicated that the new array is not assigned to the data-source object (2 is the name of the array). Can you try setting the backend to auto?

larsoner commented 7 years ago

Can you try setting the backend to auto?

Even if it fixes it, I don't think it's a good solution to the problem. This implies that these changes propagate problems to other packages, so it would be better if we could find some fix that doesn't do this.